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Prsface 


This  report  proposes  to  investigate  and  apply  various  quantal  re¬ 
sponse  models  to  determine  the  applicability  o-f  these  methods  on  data 
generated  by  the  Avionics  Evaluation  Program.  It  is  hoped  that  this 
effort  will  be  a  basis  -for  moving  quanta!  response  methods  out  of  the 
realm  of  bio-assay  and  into  more  general  applications. 

I  wish  to  acknowledge  my  indebtedness  to  my  thesis  committee,  Or. 
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ies,  for  his  support  of  this  effort,  and  particularly  for  the  computer 
source  code  and  program  deck  which  he  gave  me.  I  would  also  like  to 
thank  my  wife,  Li  Hua,  for  her  support  and  understanding  during  this 
effort. 
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Abstract 


The  United  States  Air  Force  has,  over  the  past  decade  or  so,  in¬ 
vested  much  time  and  money  in  computer  simulations  and  models.  At  the 
most  basic  level  almost  all  of  these  simulations  are  input-output  type 
procedures;  variables  of  interest  are  changed  to  determine  the  effect 
they  have  on  some  other  factor.  This  process  is  virtually  indistin¬ 
guishable  from  dose-response  problems  in  bio-assay,  hence,  is  capab’  ?  of 
being  analyized  by  the  same  methods  used  in  bio-assay.  The  two  mos 
commonly  used  techniques  are  probit  and  logit,  but  there  are  many  other 
available  techniques.  An  alternative  to  performing  numerous,  and  som».— 
times  redundant,  simulations  is  to  use  these  techniques  whenever  pos¬ 
sible. 

^ - \ 

Data  from  the  Avionics  Evaluation  Program  <AEP>  were  used  as  the 
basis  for  estimating  the  probability  of  aircraft  abort,  based  on  the 
mean-time-between-fai lure  <MTBF>  of  various  equipment  items,  using  four 
quanta!  assay  techniques.  The  fits  obtained  from  these  models  were  com¬ 
pared  to  the  more  popular  probit  and  logit  results  previously  obtained 
by  Dr.  David  Barr. 


EXPLORATION  OF  DOSE-RESPONSE  TECHNIQUES  WITH  SOME 


APPLICATIONS  TO  A  SIMULATION  PROBLEM 


I  INTRODUCTION 


Background;  An  Analysis  Problem 

It  is  common  for  analysts  to  use  computer  modeling  and  simulations 
in  problem  solving.  The  reasons  for  using  simulations  are  numerous. 
However,  simulations  usually  help  the  analyst  determine  the  response  of 
some  system  to  a  change  in  that  system's  environment,  or  operating  char¬ 
acteristics.  In  other  words,  simulations  help  the  analyst  describe,  pre¬ 
dict,  or  simply  understand  the  behavior  of  a  complex  system  under  a  given 
set  of  circumstances. 

There  is  no  doubt  that  simulations  are  useful,  but  is  it  always 
practical  or  necessary  to  use  them?  It  is  not  hard  to  think  of  situa¬ 
tions  where  the  answer  is  no,  and  the  specific  situation  described  be¬ 
low  is  just  one  such  example. 

My  last  project  as  an  analyst  for  the  Avionics  Laboratory  at  WPAFB 
was  to  run  computer  simulations  of  a  mission  analysis  program  (Appendix 
A).  The  project  was  to  determine  the  mission  effectiveness  of  the  ATF 
(advanced  tactical  fighter)  with  a  mixed  suite  of  existing  and  concep¬ 
tual  (new)  avionics.  One  of  the  measures  of  effectiveness  was  aircraft 
aborts  due  to  failure  of  a  particular  piece  of  equipment,  or  subsystem, 
with  varying  MTBFs  (mean-time-between-failure) . 


I  performed  several  thousand  simulations  for  this  analysis.  Unfor¬ 
tunately,  even  with  such  a  large  number  of  simulations  all  I  had  was  a 
dozen  or  so  data  points  (i.e.  MTBF-abort  pairs)  for  each  piece  of  equip¬ 
ment.  If  there  was  a  need  for  an  abort  rate  at  some  MTBF  not  previously 
considered  then  either  another  simulation,  or  some  estimate  based  on  the 
existing  data  would  be  required. 

Since  it  was  not,  and  usually  is  not,  practical  to  run  computer 
models  for  every  conceivable  point  it  was  obvious  that  there  was  a  need 
for  some  curve-fitting  technique.  It  was  at  this  point  that  my  working 
group  contacted  Dr.  David  Barr  (Air  Force  Institute  of  Technology)  and 
asked  him  to  study  the  problem.  He  determined  that  there  were  tech¬ 
niques  for  ‘...estimating  probabilities  when  given  a  set  of  relative 
frequencies,  each  obtained  as  the  response  of  a  system  to  a  level  of 
quantitative  stimulus,  known  as  probit  and  logit  analysis”  (Ref  3:1). 

Probit  (probability  unit)  analysis  originated  in  biology  and  its 
application  in  that  area  is  widely  accepted.  The  following  is  the  gen¬ 
eral  concept  of  probit  (parenthetical  matter  is  my  own): 

An  analyst  is  interested  in  the  effect  of  same  drug 
(failure  rate)  on  the  survival  (aborts)  of  a  large 
number  of  insects  (systems).  One  possibility  is 
that  each  insect  (system)  survives  until  a  certain 
critical  dosage  (MTBF)  is  reached  and  that  they  all 
die  (abort)  as  soon  as  this  limit  is  surpassed;  but 
that  is  an  extreme  case.  It  is  much  more  plausible 
that  the  critical  level  varies  from  one  insect  (sys¬ 
tem)  to  the  other  according  to  a  certain  distribu¬ 
tion.  Wien  there  are  many  independent  factors  de¬ 
termining  the  critical  level  for  each  insect  (sys¬ 
tem),  the  central  limit  theorem  may  be  used  to  jus¬ 
tify  the  choice  of  the  normal  distribution.  Thus, 
when  p  is  the  proportion  of  insects  (systems)  killed 
(aborted),  the  analyst  applies  the  probit  transfor¬ 
mation  y®F(p)  and  he  then  proceeds  to  express  y 
linearly  in  terms  of  the  dosage  (MTBF)  of  the  drug. 

(Ref 34  :630)  . 


It  you  were  to  plot  these  dose-response  (input-output)  curves  with 
the  dose  on  the  horizontal  axis  and  response  on  the  vertical  axis  the 
curve  is  S-shaped,  or  sigmoid,  and  is  similar  to  the  cumulative  distri¬ 
bution  function. 


...it  seems  natural  to  try  to  -fit  a  distribution 
function  to  the  points...  The  probability  curve 
which  came  to  mind  first  was  the  normal,  or  Gau¬ 
ssian,  distribution  function;  since  negative  val¬ 
ues...  have  no  physical  interpretation,  it  made 
sense  to  make  a  logarithmic  tra  ormation  to  the 
lognormal  distribution.  This  r  ulred  in  what  is 
known  as  orobi t  anal  vsis.  As  1  ar  approach, 
in  which  the  logistic  growth  c  -e  is  used  in  place 
of  the  Gaussian  distribution  f  ♦■ion,  is  known  as 
looi t  anal vsis  (Ref  3s 5)  . 

Dr.  Barr's  work  (Ref  3)  on  the  problem  of  curve-fitting  showed  po¬ 
tential  for  the  application  of  these  statistical  techniques  to  the  'in¬ 
put-output'  type  problem  discussed  earlier.  However,  while  these  tech¬ 
niques  worked  well  for  some  equipment  items,  it  worked  only  marginally 
well  for  others,  and  not  at  all  for  still  others.  This  leads  to  the 
following  question;  are  there  other  dose-response  techniques  available 
to  analyze  the  existing  data  (Appendix  C)  which  would  give  either  better 
fits,  or  at  least  fit  those  items  for  which  the  probit  and  logit  methods 
were  only  marginal?  This  question  is  the  underlying  basis  for  this  ef- 


Before  reviewing  existing  dose-response  techniques  (Section  II)  it 
is  necessary  to  state  some  constraints  of  the  existing  data  which  limit, 
or  eliminate  the  use  of  certain  techniques.  It  should  be  kept  in  mind 


that  these  constraints  are  limiting  factors  only  for  this  effort;  future 


efforts  may  be  well  advised  to  use  techniques  which  l  could  not. 

The  number  one  constraint  is  the  amount  of  existing  data,  152  data 
points  total  -for  17  equipment  items.  This  averages  out  to  about  nine 
data  points  -for  each  equipment  item,  and  -tor  some  analyses  this  would  be 
suff icient.  Unfortunately,  100'/i  o-f  this  data  yields  abort  rates  less 
then  3 X;  this  is  Known  as  the  low-dose  range  (usually  considered  <  0.10) 
and  presents  problems  o-f  its  own  -for  extrapolation  over  the  rest  o-f  the 
curve. 

While  it  would  have  been  convenient  to  run  more  simulations  to  ob¬ 
tain  a  wider  range  o-f  abort  rates,  it  is  no  longer  possible  to  do  so. 

The  AEP  model  (Appendix  A)  was  removed  from  use  just  prior  to  this  ef¬ 
fort.  But  even  before  removal  the  model  had  undergone  extensive  modifi¬ 
cation  and  enhancing,  which  would  have  made  any  comparison  of  new  and 
old  data  suspect. 

The  second  constraint,  while  not  a  problem,  eliminates  the  use  of 
techniques  known  as  mutihit  and  multistage  (Refs  2,  17:1277-1278). 

These  models  will  be  discussed  in  Section  II.  Basically,  however,  the 
problem  is  that  unlike  living  organisms,  which  can  be  exposed  to  a  sub¬ 
stance  then  periodically  re-exposed  (rehit)  until  a  tumor  or  other  re¬ 
sponse  is  obtained,  a  hardware  item,  within  the  AEP  model,  fails  based 
only  on  its  designed  (ore  time)  MTBF. 

A  third  constraint  again  stems  from  the  fact  that  performing  new 
simulations  is  not  possible,  which  eliminates  the  use  of  sequential 
methods.  As  the  name  implies,  this  technique  involves  running  a  simu¬ 
lation,  or  experiment,  and  then  performing  another  simulation  at  either 
a  higher  or  lower  level  of  the  stimulus  based  on  the  previous  results. 
You  then  repeat  this  procedure  until  obtaining  the  results  or  accuracy 


desired  (Refs  3s  1-2,  34). 

The  last  set  of  techniques  which  cannot  be  used  are  known  as  'time 
to  occurrence,  or  time  to  response*  (Refs  9:161-143,  17:1284-1236). 
Since  data  is  not  available  concerning  when  items  failed  there  is  no 
ready  data  base  on  which  to  test  these  techniques. 


»mal  Probli 


?nt  anc 


In  a  preliminary  analysis  effort  performed  by  myself  (Ref  26)  for 
the  Avionics  Laboratory  it  was  of  interest  to  determine  the  effect  of 
various  avionics  equipment  MTBFs  on  the  abort  rates  for  the  ATF.  Due  to 
time  considerations,  complexity  of  the  simulation  model,  and  working 
group  resources  it  was  desired  to  find  some  mathematical  tec;  e  for 
predicting  abort  rates  over  a  range  of  MTBFs,  using  a  limited  number  of 
simulations  for  each  equipment  item. 

A  study  conducted  by  Dr.  Barr  (Ref  3)  showed  the  applicability  of 
using  quantal  (dose-response)  assay  techniques  for  determining  these  a- 
bort  rates.  In  particular,  he  was  able  to  fit  most  of  the  17  equipment 
items  under  consideration  using  the  probit  and  logit  models. 

It  is  the  purpose  of  this  effort  to  determine  what  other  quantal 
assay  methods  exist  which  may  fit  the  avionics  equipment  in  question  to 
dose-response  curves.  More  specifically:  what  quantal  models  are  avail¬ 
able  that  will  estimate  the  probability  of  an  abort  given  an  MTBF  for  a 
particular  piece  of  avionics  equipment  using  the  existing  data  of  Appen¬ 
dix  C?  The  answers  to  this  question  have  implications  for  analysts  in 
general,  particularly  where  the  analysts7  situation  involves  the  use  of 
simulations. 


Section  II  contains  a  summary  of  the  types  of  existing  quantal 
< dose-response)  models.  A  brief  description  of  each  type  is  given  a- 
long  with  the  mathematical  development  where  applicable. 

In  Sections  III  thru  VI  the  one-hit,  a  generalization  of  the  probit 
and  logit,  quant  it,  and  a  symmetric  and  asymmetric  transformation  are 
used  to  analyze  and  fit  the  existing  sample  data.  In  these  sections  the 
models  are  explained  along  with  techniques  for  implementation.  Also 


discussed  are  special  considerations  and  limitations  of  these  models  as 
well  as  results  and  interpretation  of  the  results. 


i  •  . 


II  Review  of  Dose-Resoonse  Methods 


Any  system  that  yields  a  response  to  a  given  stimulus  can  use  the 
methods  listed  below.  The  almost  exclusive  use  o-f  these  methods  in  bio¬ 
assay  is  no  reason  for  their  non-use  in  other  areas.  The  only  real  con¬ 
straint  is  that  the  system  response  variable  must  be  a  random  variable 
that  takes  on  only  one  o-f  two  values;  success  or  failure,  0  or  1,  yes  or 
no,  etc.,  ‘...such  observations  are  called  binary;  an  older  term  is 
ouantal ‘  (Ref  12; 1) . 

There  are  numerous  quantitative  theories  that  attempt  to  relate  the 
frequency  of  response  to  the  level  of  stimulus.  Crump  (Ref  ?>  and  Fish- 
bein  (Ref  17)  categorize  the  methods  most  commonly  used  into  two  major 
types:  dichotomous  response  models  and  time-to-response  models. 

QlChotomovs_RgSBgn?e  Made) 3 

One-Hit  Model  and  Extensions.  The  most  elementary  dose-response 
model  is  the  one-hit,  or  linear,  model.  ‘The  one-hit  model  is  obtained 
by  assuming  that,  with  the  exception  of  Xd  hits  at  dosage  d,  the  proba¬ 
bility  of  exactly  x  hits  is  given  by  the  general  term  of  the  Poisson 
distribution...  *  (Refs  17:1277,  33).  The  general  term  is  as  follows 

P(X»x)  ■  (Xd)XEexp(-xd)3/x!  (2.1) 

Clearly  if,  as  in  our  case,  only  one  hit  is  required  to  produce  a  re¬ 
sponse  then  the  above  col  1  apses  to 

P(d) *P(x  > 1>  =  1-P(x< 1)  =  l-exp(-Xd)  (2.2) 


n 


•  «  *  «  •  " 


where  X  is  the  process  rate,  or  rate  of  change  of  the  dose  response 
curve  at  d=0. 

The  point  was  made  by  Fishbein  (Ref  1?)  that  if  you  are  working  in 
the  low  dose  region  (p<0.1)  then  xd  is  small  and  P(d)=xd.  This  implies 
a  simple  linear  model  where  the  response  is  directly  proportional  to  the 
dose,  with  slope  X.  Since  all  of  the  sample  data  (KTBF-abort  pairs)  is 
well  within  the  low  dose  region  this  is  one  possible  model  to  use. 

The  natural  extension  to  the  one-hit  model  is  the  multihit  model. 
This  model  considers  that  i-f  at  least  K  hits  of  a  dose  d  are  required  to 
produce  a  response  then  (Ref  17:1277) 


P(d)  =  1 


K-l 

-  V  (xd)  lexp(- 


Note,  if  K  is  allowed  to  take  on  non-integer  values  then 


P(d)«  /nytk"1exo(-t)  dt 
Jq  (k-l)! 


Rewriting  this  as: 


P<d)«P<d 


i*,x)»  f 


d.ki.K-l 


exp(-Xt)  dt 

r<k) 


yields  the  generalized  multihit  dose-response  model,  or  more  simply  a 
gamma  distribution  with  scale  parameter  X  *  id  a  shape  parameter  K  (Ref 
33:342) . 

Another  extension  or  further  generalization,  if  you  like,  of  this 
stochastic  process  is  the  multistage  model,  ”...  where  the  lifetime 
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probability  of  tumor  induction  can  bo  expressed  as 


P(d)-i-exp[-<a1-f1dX<r2-f2d)  . . .  ( c^-^d)  3  (2.3) 

where  a. 26,  *^>0,  and  K  represents  the  number  of  transitions  or  muta¬ 
tional  stages  in  the  carcinogenic  process"  (Refs  17:1278,  2).  This 
model  has  no  application  to  our  hardware  items  since  the  dose  (MTBF)  is 
applied  only  once  (as  an  input  parameter,  which  remains  fixed,  to  the 
AEP  model,  see  Appendix  A),  and  where  the  above  would  collapse  to: 

P(d)»l-exp[-(a-Pd) ] 

However,  the  multistage  mode)  could  be  useful  if  you  thought  about 
failures  of  an  item  which  did  not  cause  an  abort  as  the  mutation.  Then 
after  repairing  the  failed  item  it  is  replaced  in  the  aircraft  and  the 
aircraft  is  flown  again  (note  that  the  repair  restores  the  MTBF  and  is 
the  'rehit'  of  the  item).  This  process  is  repeated  until  a  failure  (the 
number  of  failures  would  have  to  be  determined  by  some  stochastic  pro¬ 
cess)  causes  an  abort  (tumor). 

In  an  article  by  Guess  and  Crump  (Ref  19)  there  is  one  multistage 
model  worth  separate  discussion.  It  is  a  general  polynomial  model  for 
dealing  with  low-dose  extrapolation,  and  is  the  only  multistage  model 
which  is  well  documented  and  supported  in  a  series  of  articles  (by  Guess 
and  Crump) . 

To  obtain  this  model  first  consider  the  following: 


P(d)«l-expt- 


[a>P.  (d>‘ 


...it  is  assumed  that  K  >1  different  events  must  oc¬ 
cur  in  a  single  cell  berore  a  cancer  (response)  is 


initiated  ...  a  is  the  rate  at  which  an  initiating 
event  for  the  ith  stage  occurs  due  to  a  spontaneous 
(i.e.  not  dose  related)  carcinogenesis,  and  ?.  is  the 
per  JL  th  power  dose  at  which  an  initiating  ev&nt  for 
the  ith  stage  occurs  due  to  dose-related  carcinogen¬ 
esis  (Ref  19:17). 


Now  rewriting  the  above  the  model  becomes 


P#(d)sP(d)»l-exp(-f(d>)  (2.4) 

where 

K 

f(d)*^f.d1  f.>0 

i*0 

and  where  f  is  a  polynomial  with  nonnegative  coefficients.  K  is  the  de¬ 
gree  of  f,  and  along  with  the  coefficients  of  f  must  be  estimated  (Ref 
19:18). 

Note  the  similarity  of  Eq  (2.4)  with  the  general  multistage  model 
in  Eq  (2.3).  This  model,  however,  uses  a  general  polynomial  of  unknown 
degree  to  fit  the  data.  This  model  is  of  particular  interest  if  one  is 
working  in  the  medium,  or  high-dose  range  and  then  wishes  to  extrapolate 
to  the  low-dose  range  for  risk  estimates  (Ref  19:21-22). 

To  construct  confidence  intervals  Crump  et.  al .  have 


...developed  'envelope  curves'  which  are  constructed 
for  both  risk  and  dose  ...  these  curves  are  con¬ 
structed  by  binomial ly  simulating  188  sets  of  dicho¬ 
tomous  dose-response  data,  representing  100  indepen¬ 
dent  replications  of  the  same  experiment  ...  (with) 
the  same  set  of  test  doses...  (Ref  10:440). 


Guess  and  Crump  (Ref  20)  also  develop  maximum  likelihood  estimation 
techniques  for  this  polynomial  model. 


methods  of  analysis  of  binary  data  are  probit  (Ref  15)  and  logit  (Ref  5, 
12).  For  the  probit  (log-probit)  model  the  probability  of  a  response 
being  induced  by  a  stimulus  (dose)  d  is  given  by 

P(d)  =  $(<x+  flog  d)  (2.3) 

where  0  denotes  the  standard  cumulative  Gaussian  (normal)  distribution. 

The  logit  model  like  the  probit  model  leads  to  an  S-shaped  dose- 
response  curve;  its  equation  is 

P(d)  »  1/1  l+exp(-(a+P1og  d))l  (2.6) 

It  approaches  zero  response  as  dose  d  decreases  more  slowly  than  does 
the  probit  curves  since  I im(P(d)/dP)=constant  as  <H0  (Re-f  17:1279). 

Dr.  Barr's  work  (Ref  3)  on  the  curve-fitting  problem,  using  probit 
and  logit,  yielded  the  values  in  Tables  I  and  II  for  goodness-of-f it 
based  on  chi-square  tail  probabilities  for  the  17  equipment  items  con¬ 
sidered.  These  values  will  be  the  bench  mark  against  which  all  other 
model  fits  will  be  tested  (primarily  since  the  probit  and  logit  models 
are  well  developed  and  widely  accepted). 

These  fits  for  the  probit  and  logit  are  not  really  very  different. 
This  is  somewhat  expected  since  according  to  Finney  the  logistic  and 
normal  distributions  are  ”...  scarcely  distinguishable  ...  between  re¬ 
sponse  rates  of  0.91  and  0.99...”  (Ref  16:496) . 

There  is  a  method  described  by  Chambers  and  Cox  (Ref  8)  which  may 
better  discriminate  between  the  logit  and  probit  models.  However,  it 
depends  on  having  a  few  dose  (MTBF)  levels  then  performing  a  test  which 


will  determine  the  appropriate  spacing  of  the  dose  (HTBF)  levels  -for 
this  discrimination.  Performing  more  simulations  is  impossible  at  this 
time,  hence  this  method  cannot  be  used. 

The  vast  majority  o-f  the  literature  encountered  in  this  review  was 
in  the  biological  sciences,  and  thus  there  was  a  propensity  to  find  the 
LD50.  LO50  is  Known  as  the  59’/.  lethal  dose,  or  the  median  effective 
dose.  While  this  effort  has  no  interest  in  the  LD50  <LD01  or  less  would 
be  more  informative  for  design  engineers)  it  is  worth  mentioning  since 
it  brings  up  the  matter  of  transformations. 

Finding  the  estimates  of  F,  and  02  for  the  distribution  given  by 
Eq<2.3)  is  generally  by  means  of  the  probi t  transformation  of  the  exper¬ 
imental  results.  "The  probit  of  the  proportion  P  is  defined  as  the  ab¬ 
scissa  which  corresponds  to  a  probability  P  in  a  normal  distribution 
with  mean  3  and  variance  1*  (Ref  13:21).  That  is,  the  probit  of  P  is  Y, 
where 

-1/2  /*Y’5 

P-  <210  */  exp-(u*/2)  du  (2.7) 

This  transformation  from  proportions  to  probits  has  the  effect  of 
straightening  out  the  norma)  S-shaped  curve.  Comparing  Eqs  (2.5)  and 
(2.7)  shows  that  the  probit  Y  is  related  to  dosage  d  by  the  simple  lin¬ 
ear  equation  Y«5+(d-l»)/d  (Ref  15),  and  now  to  estimate  LD50  you  simply 
find  the  value  of  d  which  gives  Y"5.  The  usefulness  of  this  transforma¬ 
tion  is  to  simplify  mathematical  calculations. 

There  are  numerous  transformation  techniques  available  (Refs  4,  6, 
22).  They  depend  solely  on  the  models  used,  the  experimental  data,  comp¬ 
utational  considerations,  convenience  of  their  use,  and  ability  to  in- 


terpret  results  meaningful  1 y. 


A  generalization  of  the  probit  and  logit  models  is  mentioned  in  an 
earlier  paper  by  Prentice  (Ref  31)  and  completely  detailed  by  Prentice 
in  a  latter  paper  (Ref  30).  This  model  takes  the  form  of 


P(d) 


■/ 

•r—CC 


f  (w)  dw 


(2.8) 


where  /“(d-IO/d  and  are  to  be  estimated.  The  pdf,  f (w) ,  has  the 
following  form 


f(w)  =  exp(wm) ( 1+expw) (  m  n)  (2.9) 

£(m,n) 

where  £  is  the  beta  function.  The  logistic  model,  Eq  (2.6),  is  given  by 
m*n=l  and  converges  to  the  normal  distribution  as  m,n-*w.  Other  special 
cases  for  various  m,n  values  are  also  given  (Ref  30:762).  This  model 
will  be  discussed  more  fully  in  Section  IV  where  it  is  applied  to  the 
sample  data. 

Note  this  model  is  really  nothing  more  then  a  beta  of  the  second 
kind  distribution  with  the  transformation  u*exp(w> ,  0<u.  To  show  this 
note  that 


f (u)«f (w) Idw/du | 

■  u*/£  ( l+u)  *m+n*J(m,n)  ul 

(m—l).,/ti  .  (m+n) „  > , 

*  u  /£(l*u)  *(m,n)3 


or  a  beta  of  the  second  kind 


Quant  i  t  Analysis.  In  a  paper  by  Copenhaver  and  Mi  Ike  (Ref  11)  a 
new  technique  is  offered  for  analyzing  quantal  responses,  which  they 


call  Quantit  analysis.  The  underlying  distribution  is  called  the  omega 
distribution  by  Copenhaver  and  Mi  Ike  and  is  characterized  by  a  cdf  of 

F  ( x  ( q) )  =q  (2.10) 


and  a  pdf  of 


f(x(q))=l-|2q-l|v+1 


where  0<q<l  and  v>-l  ,  and  where 


x 


<q> 


-L 


dz/f <x(z) > 
1/2 


(2.11) 


As  noted  by  the  authors  this  distribution  is  a  double  exponential  when 
v*0,  logistic  when  v-1,  and  uniform  for  the  limiting  case  as  v-»*. 

If  we  let  P  be  the  probability  (or  proportion)  of  response(s)  at 
dosage  x.f  and  we  let  P.=F(a+Px.)  then  the  tolerance  distribution  is 
given  by 


and 


f  (cr+£x  .  )  =  l-|2p. -1 

l  l 


.v+1 


(2.12) 


a+Px.*h  (p.)*/* 'l/(  l-|2z-l 
i  v  i  ^,/2 


v+1 

I  )  dz 


(2.13) 


hv(pi>  is  termed  the  'quantit*  of  p. (Ref  11:178).  The  computations  to 
obtain  the  parameters  are  similar  to  that  for  the  model  presented  by 
Prentice  which  is  discussed  in  Section  IV.  Section  V,  however,  will 
discuss  more  fully  the  individual  computations  for  the  Quantit  model, 


and  the  application  of  it  to  our  sample  data. 


Aranda-Ordaz  Family  of  Trans-forms.  One  of  the  reasons  for  first 
using  the  logit  model  in  the  analysis  of  binary  data  is  its  simplicity, 
since  the  logistic  transformation  used  to  obtain  the  logit  method  is  a 
linear  function  of  the  parameters. 

The  logistic  transform  is  simple,  but  is  it  adequate?  Aranda-Ordaz 
(Ref  1)  deve lopes  simple  test  procedures  to  determine  both  symmetric  and 
asymmetric  departures  from  the  logistic  distribution.  He  also  develops 
two  families  of  power  transforms  which  are  alternatives  to  the  logistic. 
These  symmetric  and  asymmetric  alternatives  each  include  the  logistic  as 
a  special  cases. 

The  symmeteric  family,  which  essentially  yields  the  same  results 
when  working  with  either  successes  or  failures  (Ref  1:358),  is  given  by 
Equation  (2.14) 

T  (9)  -  (2/X)  ( 9X-(1-6)X)  (2.14) 

\ex*(i-8>y 

and  in  the  limit  is  the  logistic  when  X*0  and  a  simple  linear  for  x=l. 

Now  solving  Eq  (2.14)  as  a  function  of  T  yields: 


9<D- 


’  e 

_ ( i+vr/2) 1/x 

I  (  l*XT/2)  1/x+  ( 1-X1/2)  1/x 


(*Xf  <  -1) 

1X1/2 K  1  (2.15) 

(#X1  >  1) 


We  assume  that  T  has  a  1  inear  expression  in  terms  of 
some  parameters  associated  with  the  explanatory  var¬ 
iables  considered  in  a  specific  situation...  If  we 


VN 


•fit  by  maximum  likelihood  a  linear  expression  for  t 
for  a  range  of  values  of  X  we  may  consider  the  max¬ 
imized  log  likelihood  as  a  function  of  x  and  hence 
derive  not  only  the  maximum  likelihood  estimate 
but  also  determine  which  values  of  provide  an  ac¬ 
ceptable  fit.  (Ref  1:358). 

The  aysmmetric  model  is  given  by  the  family  (assuming  log  W(8)=T) : 

Wx(e>  »  ((l-8)”X-l)/x  (2.16) 

This  is  the  logistic  model  when  X»l,  and  again  solving  as  a  function  of 
T  we  get 

f 

l-(  1+X  exp(T))  1/X  (X  exp(T))>-l 

e(T>»  . 

1  otherwise 

The  same  assumption  about  T,  and  the  same  procedure  as  above  is  used  to 
obtain  values  for  X. 

There  are  tests  (Ref  1:368-361)  to  determine  if  there  are  any  sym¬ 
metric  or  asymmetric  departures  from  the  logistic.  The  attractiveness 
of  these  tests  are  that  they  can  be  conducted  using  "...(the  values) 
computed  from  the  output  of  the  logistics  fit.*  (Ref  1:360).  Another 
by-product  of  these  tests,  specifically  the  symmetric  test,  is  that  it 
mav  permit  "...  discrimination  between  the  logistic  and  probit  mod¬ 
els...*  (Ref  1:361). 

The  above  models  and  tests  will  be  fully  discussed  and  implement ted 
in  Section  VI.  In  that  section  I  will  try  to  discriminate  between  the 
logistic  and  probit  fits  already  performed  by  Dr.  Barr  (Ref  3). 


Time  to  Occurrence 


According  to  Fishbein  'The  second  type  of  dose-response  modeling 
that  is  receiving  increasing  attention  deals  with  the  distribution  of 
the  'time  to  occurrence'  (latent  period)  and  its  relation  to  dose.'  (Ref 
17:1284).  Unfortunately  the  few  models  which  were  discussed  by  Fishbein 
(Ref  17)  and  Crump  (Ref  9)  all  had  criticisms  leveled  at  them  by  the  au¬ 
thors  and  others.  Since  these  models  are  not  well  accepted,  and  since 
(as  in  the  multistage  models)  new  interpretations  of  the  dose-response 
processes  are  needed  this  type  is  mentioned  only  for  completeness. 

In  general  this  type  of  model  could  be  defined  as  'the  time  of 
death  (abort)  from  the  type  of  cancer  (failure)  of  interest  or  as  the 
time  of  the  first  appearance  or  detection  of  a  particular  tumor  type' 
(Ref  17:1284),  parenthetical  matter  my  own. 


Ill  The  One-Hit  Model 

The  simplest  curve  -fitting  technique,  i-f  we  exclude  drawing  a  line 
between  two  points,  is  the  simple  linear  regression  model.  The  one-hit 
model  is  nothing  more  than  a  linear  regression  model  with  a  logarithmic 
transform  of  the  data.  Regression  is  a  well  understood  and  straightfor¬ 
ward  technique,  and  as  such  needs  no  separate  discussion,  except  for  one 
special  consideration  discussed  later  (regression  through  the  origin). 

If  we  assume  that  it  takes  exactly  x  failures  (hits)  of  some  item 
to  cause  an  abort  for  a  particular  MTBF  (dose),  where  the  failures  are 
independent,  random  events,  then  we  can  use  the  Possion  process  Eq  (2.1) 
to  describe  the  probability  of  an  abort.  However,  as  mentioned  before, 
if  it  takes  only  one  failure  to  cause  an  abort  then  the  general  Possion 
term  collapses  to  Eq  (2.2),  or  P(d>  =  l-exp(-xd)  where  X  is  the  pro¬ 
cess  rate. 

It  is  easy  to  verify  that  if  P  is  small  (P  <  0.1)  then  P(d)*xd. 
This  implies  that  for  small  P  a  linear  equation  going  through  the  origin 
describes  the  data.  While  our  data  is  much  less  then  10X,  it  comes  from 
high  MTBFs  (doses).  This  situation  is  exactly  the  opposite  of  the  usual 
bio-assay  problem  where  a  low  dose  causes  a  low  response  rate. 

Looking  at  Figure  1  shows  that  it  makes  no  sense  to  force  the  known 
portion  of  the  curve  (solid  line  segment)  through  the  origin.  This  is 
because  this  portion  of  the  curve  has  a  negative  slope.  If  we  flip  the 
curve  around,  as  shown  in  Figure  2,  then  we  might  be  able  to  force  the 
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known  portion  of  the  curve  through  the  origin  and  use  the  simple  equa¬ 
tion  P(d)=Xd.  However,  we  cannot  use  this  method  since  P  is  no  longer 
small  and,  hence,  not  approximated  by  xd.  This  is  no  handicap  though, 
since  Eq  (2.2)  is  not  hard  to  apply  directly. 

Now  if  we  let  P  be  the  proportion  of  non-aborts  to  launches  then 
Q*l-P,  and  Eq  (2.2)  can  be  written  as: 

Q  *  exp(-xd)  (3. 1) 

and  taking  logs  of  both  sides  yields  Lnq  =  -Xd,  or 

-Lnq  =  Xd  (3.2) 

To  apply  this  linear  model  means  forcing  the  equation  through  the 
origin,  since  there  is  no  constant  term.  The  Control  Data  Corporation 
Cyber  750  computer  implementation  of  the  Statistical  Package  for  the 
Social  Sciences  (SPSS)  (SPSS  is  widely  available)  allows  one  to  force 
an  equation  through  the  origin  using  the  appropriate  option  (option  1?) 
in  the  regression  procedure.  However,  each  of  the  correlation  coeffi¬ 
cients,  R2  and  adjusted  R2,  are  unadjusted  for  the  mean  when  using  op¬ 
tion  19.  But,  SPSS  displays  an  extra  line  of  output  with  these  values 
adjusted  for  the  mean  as  suggested  by  Theil  (Ref  34.. :6) . 

When  forcing  an  equation  through  the  origin  using  SPSS  one  should 
first  determine  the  appropriateness  of  this  option.  In  our  case  it 
makes  perfect  sense,  since  we  would  expect  to  have  a  'continuous-  abort' 
for  an  MTBF  of  zero  for  any  item  which,  by  itself,  causes  an  abort.  An¬ 
other  item  for  consideration  is  the  adjusted  for  the  mean  correlation 


coefficients,  R*  and  adjusted  R*.  If  these  values  should  happen  to  go 
negative  then  forcing  the  equation  through  the  origin  is  not  appropriate 
(Ref  34 : 177). 

After  performing  an  initial  regression  on  all  17  equipment  items 
(Appendix  C> ,  and  observing  that  the  correlation  coefficients  adjusted 
for  the  mean  were  al I  negative,  a  new  approach  was  undertaken.  In¬ 
stead  of  using  the  MTBF  the  simple  transform  d'»Ln  MTBF  was  tried,  or 
-Ln<f«xd'»XLnd.  If  we  now  solve  for  P  we  obtain: 

P  «  l-(d)"X  (3.3) 

With  Eq  (3.3)  note  that  the  value  of  the  MTBF  can  approach  but  not 
equal  zero.  Since  there  is  no  reason  not  to  allow  the  MTBF  to  equal 
zero,  at  least  in  theory  if  not  actual  practice,  I  applied  another  sim¬ 
ple  transformation,  really  a  translation.  The  translation  was  to  simply 
add  one  to  the  MTBF  before  taking  logs,  yielding: 

P  ■  l-exp(-\  Ln(d«-1>> 
or 

P  -  1-(CH1)’X  (3.4) 

Besides  letting  the  MTBF  take  on  all  non-negative  values  Eq 
(3.4)  has  the  property  of  being  a  known  distribution,  the  Pareto  distri¬ 
bution,  translated  by  a  value  of  one. 

Using  Eq  (3.4)  and  performing  regression  on  all  17  equipment  items 
gave  acceptable  fits  for  all  of  the  items.  The  determination  of  accept¬ 
ability  was  by  the  usual  methods;  checking  the  F-test  values,  correla- 


tion  corf f icients,  residuals,  etc.  Listed  below  are  the  equations  ob¬ 
tained  from  the  regression  procedure  (remember  P  is  the  proportion  of 
non-aborts) . 
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P  » 

1-expt -.9442832?  Ln<d*l>) 
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P  ■ 

i-expt- 1.9278328 

Ln(d+ 1) 1 

CDRT 

P  ■ 
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P  ■ 

1-expt -.84847599 
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P  - 
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P  » 

1-expt -.84338 133 

Ln(d*  1)  3 

SMRT 

P  » 

1-expt -.9742 1848 

Ln(d+  1) 1 

Note  the  clear  trend  in  the  above  equations  for  the  coefficient 
to  cluster  about  the  value  of  x«l.  Dr  Barr  (Ref  3)  using  probit  and  lo¬ 
git  had  no  such  trend  for  any  of  his  coefficients. 


in  of  the  Fit  and  Normality  Assumptions 


Determinatio 


During  the  preliminary  AEP  analysis  (Ref  26)  there  was  a  real  con¬ 
cern  about  the  random  number  generator  used  in  the  AEP  model.  For  a  -few 
o-f  the  equipment  items  in  Appendix  C  there  appears  to  be  significant  de¬ 
partures  from  the  general  trend.  There  is  a  way  to  determine  if  there 
are  departures,  at  least  for  the  normality  assumptions,  when  using  re¬ 
gression.  This  is  by  performing  residual  analysis,  which  is  well  Known 
and  easy  to  perform  using  SPSS.  However,  for  a  detailed  description  of 
residuals  and  residual  analysis,  the  reader  should  see  Theil  (Ref  34). 

The  initial  SPSS  residual  plots  for  the  17  equipment  items  indi¬ 
cated  that  a  few  data  points  were  possible  'outliers'.  But,  further  in¬ 
vestigation  using  a  t<n-2)  distribution,  since  n  was  small  in  all  cases, 
revealed  no  real  outliers.  The  residual  plots  appeared  to  show  no  het- 
eroscedasticity,  or  more  simply  the  variance  appeared  to  be  constant.  I 
emphasize  the  word  appear,  since  the  residual  plots  had  as  few  as  six 
points,  and  with  so  few  points  it  would  be  misleading  to  state  that 
there  was  absolutely  no  heteroscedasticty.  However,  in  general,  there 
were  no  significant  departures  from  normality. 

To  determine  how  well  this  model  compares  with  the  probit  and  logit 
methods  I  used  the  contribution  to  the  chi-square  tail  probabilities,  as 
did  Dr.  Barr. 

If  the  true  probability  of  an  abort  for  a  given  lev¬ 
el  of  MTBF  is  Known  to  be  p,  and  if  the  number  of 
launches  is  n,  then  the  number  of  aborts  has  a  bi¬ 
nomial  distribution  with  parameter  (n,p>;  if  n  is 
large  then  the  number  of  aborts  can  be  approximated 
by  a  normal  distribution  with  mean  np  and  variance 
np(l-p).  (Ref  3:8,  also  see  14:229-230). 


Hence,  if  the  normal  approximation  is  good  for  the  abort  pro¬ 
cess  then  the  square  of  this  will  be  chi-square  distributed.  That  is, 

1/2 

if  we  let  x  be  the  number  of  aborts  then* . . .(x-np)/Cnp< 1-p) 3  will  have 

2 

a  normal  distribution...,  and  so  (x-np)  /np(l-p)  will  have  a  chi-square 

distribution..."  (Ref  3:9).  Note,  if  we  have  K  levels  of  the  MTBF  the 
chi-square  distribution  will  have  K-l  degrees  of  freedom  (df>.  We  lose 
one  degree  of  freedom  since  we  must  estimate  X. 

The  above  chi-square  distribution,  and  hence  the  validity  of  using 
chi-square  comparisons,  hinges  on  whether  or  not  the  normal  approxima¬ 
tion  of  the  abort  process  is  good.  I  bring  up  this  issue  since  the 
average  p  for  our  data  is  approximately  .002,  a  very  small  value,  and 
this  implies  the  underlying  binomial  process  is  extremely  skewed  to  the 
right.  Also  the  theory  is  for  large  n,  but  how  large  should  this  n  be? 
Our  n  was,  on  average,  approximately  5300;  is  this  large  enough? 

There  is  an  excellent  paper  by  Raff  (Ref  32)  with  some  easy  to  use 
graphs  for  determining  the  appropriateness  of  using  a  normal,  Possion, 
etc.  approximation  to  a  binomial  process.  Unfortunately  his  graphs  were 
of  little  help  in  our  case  since  p  was  so  small,  and  n  extended  beyond 
the  range  of  his  graphs. 

I  had  to  insure  that  the  normal  approximation  of  the  underlying  bi¬ 
nomial  process  was  good,  or  the  chi-square  comparisons  would  be  meaning¬ 
less.  However,  the  test  for  normality  was  simple.  First,  I  generated 
100  binomial  random  variables  with  pw.0024  and  n*5300.  Next  I  let 


y«(x-np>/tnp( 1— p> J  ,  where  x  is  the  binomial  random  variate  generated. 


Then,  using  the  Kolmogorov-Smirnov  test  in  SPSS,  I  tested  whether  y  was 
standard  normal  distributed.  The  results  were  that  I  could  find  no  evi¬ 
dence  that  the  abort  process  was  not  normal;  hence,  I  could  assume  the 
chi-square  comparisons  would  be  valid.  The  results  of  the  one-hit  chi- 
square  fit  values  is  given  in  Table  III. 


Table  III 


Chi-Sauare  Tail 

Probabi 1 i ties 

for  the  One- 

-Hit  Model 

I  tern 

Probabi 1 i tv 

Bus 

.9584 

CDRT 

.9543 

MFK 

.6789 

SLU 

.5852 

IMFK 

.4935 

MMP 

.4438 

DEK 

.3855 

MPDG 

.2127 

ART 

.1904 

SMRT 

.1313 

Processors 

.1270 

MPDS 

.0450 

INS 

.0187 

HUD 

.0095 

MTU 

.0020 

DSMU 

.0007 

SCU 

.0002 

Comparing  Tables  I  and  II  with  Table  III  it  is  seen  that  the  one- 
hit  model  gives  a  better  fit  for  the  SLU  than  either  probit  or  logit. 
Also  the  one-hit  fits  the  CDRT  and  DEK  better  than  the  probit.  However, 
with  the  exception  of  the  SLU,  no  item  in  Table  III  had  a  better  fit 
than  that  given  by  logit.  Overall,  fits  using  the  one-hit  model  were 


much  worse  than  those  obtained  by  either  logit  or  probit. 


Conclusions  and  Remarks 

The  one-hit  method  is  simple  and  easy  to  use.  Any  analyst  with 
only  a  programmable  calculator  <and  even  many  cheaper  models)  can  per¬ 
form  linear  regression.  The  methods  and  theory  of  regression  can  be 
found  in  almost  any  elementary  statistics  text,  as  well  as  many  other 
places.  However,  the  simplicity  of  the  one-hit  linear  model  should  not 
be  the  only  reason  for  its  use.  I  found  that  while  the  one-hit  model 
fit  the  data  well,  in  the  regression  sense,  it  did  not  fit  the  data  well 
by  the  criterion  of  the  chi-square  tail  probabilities.  More  simply,  at 
least  two  other  models  fit  the  data  better  than  the  one-hit. 

The  one-hit  model  is  a  valid  technique,  since  it  is  nothing  more 
than  linear  regression,  but  it  is  not  the  only  technique.  Any  analyst 
wishing  to  use  the  one-hit  model  should  insure  that  it  is  the  best  tech¬ 
nique  for  his/her  situation. 


IV  A  Generalization  of  Probit  and  Looit 


Previous  sections  indicate  that  the  most  widely  used  and  explored 
dichotomous  response  models  are  the  probit  and  logit.  In  the  literature 
reviewed  scarcely  an  article  ended  without  some  mention  of  at  least  one 
of  these  models.  However,  discrimination  between  these  two  models  is 
often  diff icul t . 


Some  consideration  (...Chambers  and  Cox  C1767]...) 
has  been  given  to  the  choice  between  probit  and  lo¬ 
git  models  with  the  general  result  that  extremely 
large  sample  sizes  are  required  to  effectively  dis¬ 
criminate  between  the  two.  Little  success,  however, 
seems  to  have  been  achieved  in  the  development  of 
sensitive  tests  of  fit  for  probit  and  logit  models 
or  in  the  development  of  alternate  classes  of  models 
when  the  usual  models  prove  inadequate  (Ref  30:761). 


Hence,  the  almost  indistinguishable  results  of  the  chi-square  fits  for 
probit  and  logit  obtained  by  Or.  Barr  (Ref  3)  are  not  surprising  (Refs 
8,  16:406) .  "This  raises  the  question  as  to  whether  tests  based  on  more 
specific  alternatives  may  be  more  sensitive*  (Ref  30:762). 

One  method  of  a  goodness  of  fit  procedure  is  to  embed  the  models 
into  a  more  general  parameteric  family  of  models  and  then  test  the  spe¬ 
cific  models,  using  ordinary  likelihood  procedures,  relative  to  the  gen¬ 
eral  one.  (Refs  7,  30).  This  is  exactly  what  Prentice  (Refs  30,  31) 
does. 

T|»f  Ifttitl 

The  probability  of  response  (abort)  for  a  given  dose  d  is  given  as 


P(d) 


■£ 


Y 

f (w)  dw 


(4.1) 
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where  y^td-IO/d  and  l*,  C  are  unknown  and  where  the  -family  is  given  by 


f(w)  »  exp(wm)  yj  exp  w) 
P(m,n) 


(4.2) 


where  P  is  the  beta  function,  and  m,n  >  0. 

2 

If  nw«n«»l  then  Eq<4.2)  is  f (w)*exp(w)/< 1  ♦  exp  w>  or  rewriting 
the  response  P(d)  *  exp(y)/(l  ♦  exp  y>  ,  which  is  a  logarithmic  trans¬ 
formation  of  the  beta  of  the  second  kind.  Prentice  (Ref  31),  after  some 
reparameterization,  shows  that  Eq(4.2)  converges  to  the  normal  distribu¬ 
tion  as  m,n-*». 


Other  limiting  special  cases  are  the  extreme  minimum 
value  <mml,n**0  and  extreme  maximum  value  (m-*w,n»l) 
densities  ...  other  limiting  distributions  are  dou¬ 
ble  exponential  <nH0,n**0>,  exponential  (m*8,n-4i>, 
and  reflected  exponential  (nH0,n*0>  (Ref  38:762). 


The  density  function  given  by  Eq(4.2)  is  symmetric  along  m*n  or 
P(d>  »  l-P(d).  Equation  (4.2)  is  also  'negatively  skewed  for  m<n  and 
positively  skewed  for  m>n"  (Ref  30:762).  Another  characteristic  of  this 
density  is  that  it  either  has  narrower  or  fatter  tails  than  the  logistic 
depending  on  whether  m>l  or  m<l  respectively  (Ref  7:1089).  'Note  that 
CEq  4.21  allows  the  choice  between  alternative  models,  such  as  the  pro¬ 
bit  and  logit  models,  to  be  reduced  to  the  choice  between  values  in  a 
single  model"  (Ref  30:762). 


As  was  the  case  for  the  one-hit,  and  quantal  response  models  in 
general,  the  fitting  of  the  distribution  to  the  resulting  sigmoid  re 


sponse  curve  is  based  on  the  conditional  probability  of  the  number  of 
responses  x  given  dose  d,  or: 


P(x  Id)  »  Q  P5*  <  1-P>n' 


(4.3) 


where  d  is  the  dosage,  n  is  the  number  of  individuals  at  dose  d,  x  is 
the  number  of  positive  responses,  and  P  is  given  by  the  distribution 
function. 

The  likelihood  function  of  Eq(4.3)  is 


k  /n.\  x .  n .-x . 

«/7  P  1  <1-P>  1  1 


(4.4) 


Now  suppose  we  have  x.  responses  for  n.  individuals  at  dose  d^ 
i»l,...,k.  Also  let  P(d)  represent  the  probability  of  response  at 
dose  dj,  where  P  depends  on  0=  ($^,...,0^),  then  the  log-likelihood 
for  0  is  simply 

k 

J.  3  log  P(d.)  +  (n^x.)  logGKd.)] 

isl 

where  GM-P. 


The  derivative  s  3  d|/d0  has  jth  component 


dl/d0.  ■  5Ztx./P(d. )  -  (n.-x.)/Q(d.)]  [dP(d.)/d0] 
J  i  i  ii  i  i  J 


(4.5) 


and  the  <j,h)  element  of  the  information  matrix  is 


T'  In./P(d.)Q(d. ) ]  IdP(d.)/de.]  tdP(d.)/d0.3  (4 .6) 

r-i  ill  ij  ih 


The  information  (Fisher  information)  matrix  is  the  matrix  of  negative 


expected  values  of  all  second  partial s  of  the  log-likelihood  function  or 


I  n<6)  =  -E^l'fxie)] 

For  more  about  the  information  matrix  the  reader  is  referred  to  DeGroot 
(Ref  14)  section  7.8,  and  Theil  (Ref  34)  section  8.4. 

As  noted  by  Prentice  calculations  for  and  its  derivatives  are 
simple  provided  P(d)  and  dP(d)/dd  are  easy  to  compute.  If  P(d)  is  given 
by  Eq(4.1)  then  a  convenient  method  for  computing  P(d),  with  underlying 
density  given  by  Eq(4.2) ,  is  1(2;  m,n)/£(m,n>  where  z=exp(y>/( 1+exp  y) 
and  I  represents  the  incomplete  beta  integral  (Ref  30:763).  If  we  let 
j0^(P ,d,m,n)  then  *  the  derivative  of  P(d)  with  respect  to  0  has  straight¬ 
forward  components  dP(d)/dF  “0  *f(y>  and  dP(d)/dd  =  <J  (Ref  30: 

763) .  However : 


and  also 


dP(d)/dm  = 


y 

(d  logf(w)/dm)  f(w>  dw 


dP(d)/dn 


■/ 

•Cod 


(d  logf(w)/dn>  f(w>  dw 


(4.7) 


(4.3) 


may  not  have  closed  form  solutions.  Prentice  suggests  that  'for  any 
fixed  (m,n)  a  straightforward  Newton-Raphson  procedure  can  be  used  to 
compute  Fr4*(m,n)  and  fe*$(m,n)'  (Ref  30:763). 

Performing  the  calculations  for  1*  and  d  as  functions  of  m,n  are,  at 
least  in  theory,  quite  simple.  First  consider 


O'*  (0,,0O)  where  Of®  (F,d)  and  0£*  (m,n> 


and 


L\  ' 


\ 

% 

B 

£ 


s/=  <Sj  ,s2>  where  s'=  (df/dl>  ,dj|/dtf>  and  s!=  (dX/dm,d(^/dn) 

Next  partition  the  in-formation  matrix  1(6)  so  that  1^(6)  is  the  upper 
2  by  2  block  o-f  1(6).  Next  -for  f ixed  6°  “a  Newton-Raphson  procedure 
iteratively  updates  trial  value  ,6°  to  £0  -H^O)-1^,  where  1^(6)"* 
and  ^  are  evaluated  at  60  ,  6°  until  convergence  to  the  MLE  is 
reached*  (Re-f  39:744).  Note  that  the  asymptotic  covariance  matrix  of 
6,  is  I  ..(6)  evaluated  at  (§.,$£)• 

We  now  have  a  method  for  determining,  by  MLE's,  the  parameters 
(i*,0).  However,  as  noted  before,  difficulty  in  calculating  P(d)  with 
respect  to  either  m  or  n  hampers  "the  use  of  asymptotic  likelihood 
methods  for  simultaneous  inference  on  all  four  parameters  (l*,d,m,n)“ 

(Ref  39:  746).  The  use  of  a  grid  of  (m,n>  values  may  overcome  this 
difficulty.  That  is,  you  maximize  the  log-likelihood  for  numerous 
fixed  values  of  (m,n) . 

There  are  same  convenient  three  parameter  submodels  of  Eq(4.2>  for 
which  the  estimation  is  more  direct.  The  specific  case  of  n=l  yields 

P(d)  =  texp(y)/(l  +  exp  y) ]m  (4.9) 

the  derivative  of  which  is 

dP(d)/dm  =  P(d)  logCexp(y)/(  1  +  exp  y)  ] 

Now  simultaneous  inference  on  all  three  parameters  (P,d,m)  is  easily 
performed  by  a  slight  modification  to  the  previously  mentioned  iterative 
procedure. 


* 


ft 
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In  theory  any  of  the  above  procedures,  the  grid  search  over  (m,n), 
the  three  parameter  model,  or  even  simultaneous  solutions  for  all  -four 
parameters  are  possible  using  a  Newton-Raphson  procedure.  However,  the 
Newton  procedure,  in  general,  has  problems  of  its  own,  *  if  the  initial 
starting  point  is  /too  far'  from  a  minimum,  the  method  will  not  converge" 
(Ref  27:443).  The  problem  of  initial  estimates  is  widely  Known,  but  in 
our  case,  as  I  will  explain,  it  is  very  evident. 

In  the  article  by  Prentice  (Ref  38)  he  applied  the  three  parameter 
model  given  by  Eq(4.9)  to  some  classical  insect  mortality  data.  The 
model  solutions  for  P,d,  and  m  were  given  and  this  gave  me  an  excellent 
chance  to  verify  and  validate  my  computerization  of  the  model  in  £q(4.9) . 
However,  I  soon  discovered  how  sensitive  to  initial  starting  values  this 
method  really  is. 

If  I  fixed  any  two  of  three  parameters  (H,d, m)  to  the  results  listed 
by  Prentice  and  then  varied  the  third,  I  found  that  I  could  only  change 
this  parameter  by  something  less  than  18X.  Outside  of  this  18%  range  the 
model  'blew-up'.  By  this  I  mean  that  the  computer  (  a  Control  Data  Cor¬ 
poration  Cyber  758)  would  either  go  into  machine  underflow  or  overflow, 
or  when  trying  to  invert  the  information  matrix  I  found  that  the  matrix 
was  ill-conditioned.  These  were  just  a  few  of  the  problems.  Remember 
this  was  just  for  one  parameter;  the  problem  intensified  directly  with 
the  number  of  parameters  varied.  This  presented  a  real  problem  since  I 
had  no  idea  as  to  the  values  of  I1  or  d  for  any  given  m  within  our  sample 
data,  let  alone  the  value  for  m. 

My  initial  reaction  was  to  do  a  search  for  m,  incrementing  it  a 


little  at  a  time;  this  would  only  leave  the  problem  o-f  estimating  l*  and 
d.  But,  as  I  soon  discovered,  different  values  o-f  the  shape  parameter  m 
required  different  initial  estimates  'for  y  and  d.  It  did  not  take  long 
for  me  to  give  up  on  this  approach. 

The  next  approach  was  the  most  costly  with  respect  to  time.  The 
approach  was  to  systematical ly  patch  the  computer  code  for  each  problem 
encountered.  I  inserted  numerous  error  checks,  put  in  bounds  to  prevent 
computer  over/underflow,  I  even  modified  the  technique  such  that  if  the 
search  procedure  started  to  diverge,  it  would  automatically  move  back 
half  the  distance  (towards  the  last  good  solution).  Unfortunately,  all 
of  this  simply  slowed  down  the  divergence,  but  the  results  were  the 
same. 

In  hindsight  the  solution  was  simple.  Since  the  problem  was  one  of 
finding  the  zeros  of  the  derivatives  of  the  log-likelihood  function,  I 
only  needed  a  routine  that  solved  nonlinear  equations.  The  solutions 
from  this  routine  would  then  be  given  to  the  Newton  procedure  as  initial 
estimates.  After  examining  many  such  techniques  I  found  one  in  the  IMSL 
package  (2XC6R)  which  worked  well.  It  is  a  conjugate  gradient  algorithm 
for  finding  the  minimum  of  a  function. 

The  conjugate  gradient  technique  is  much  more  forgiving  for  'bad' 
initial  quesses,  but  it  too  has  limits  for  these  guesses.  So  one  more 
technique  was  added  to  the  chain  of  solutions.  Note,  if  y'=a+b  In  x  ® 
(In  x-H>/d  then  l/b»d  and  M»-da,  but  y'=a+b  In  x  is  simply  the  linear 
regression  model  for  the  logit  procedure.  Using  the  regression  coef¬ 
ficients  from  the  logit  procedure,  to  obtain  initial  estimates  for  P,d, 
for  the  gradient  technique  gave  the  best  results  thus  far. 


Results  for  the  Three  Parameter  Model 

The  first  step  was  to  run  a  simple  linear  regression  on  the  log- 
linear  equation  for  logit.  Then,  using  the  coefficients  from  the  re¬ 
gression,  calculation  of  initial  estimates  for  M,d  was  performed  and 

given  to  the  gradient  algorithm  along  with  an  estimate  of  itifI  .  The 
A  A  A 

results  for  l*,d,m  out  of  the  gradient  procedure  were  then  given  to  the 
Newton  procedure.  (See  Appendix  D  for  the  computer  listings  for  the 
gradient  and  Newton  procedures.)  The  results  of  the  analysis  follow. 

In  six  cases  (ART,  INS,  MFK,  MMP,  and  SMRT)  the  initial  estimates 
out  of  the  gradient  search  algorithm  let  the  Newton  procedure  converge. 

A  A  A 

However,  the  estimates  for  F,d,m  out  of  the  Newton  procedure  were  vir¬ 
tually  unchanged  from  those  input  by  the  gradient  procedure. 

The  equation  for  the  three  parameter  model  is  given  by  Eq(4.9), 

aaa 

where  y-(d-P)/d.  The  values  of  P,d,m  for  the  six  equipment  items  which 
converged  are  given  below,  and  the  chi-square  values  are  given  in  Table 
IV. 


ART 

1**4.93807 

d=».  33509 

IDF. 01631 

DSMU 

M-6.31593 

d*. 589 15 

mF. 00487 

INS 

1**5.20106 

<t*.  81942 

mF.0 1643 

MFK 

1**5.75742 

d*. 68732 

(TF.  00448 

MMP 

1**6.83002 

d». 72084 

mF. 00344 

SMRT 


1**3.5372 


d". 90229 


mF. 85481 


Table  IV 


Chi-Souare  Tail  Probabilities  -for  the 
Convergent  Three  Parameter  Model 


Item 

MFK 

MMP 

SMRT 

ART 

DSMU 

INS 


Ml 

.8364 

.7425 

.7239 

.3871 

.2786 

.8656 


For  the  other  eleven  cases  the  estimates  out  of  the  gradient  pro¬ 
cedure  caused  the  Newton  procedure  to  diverge.  Even  though  the  esti- 

A  A  A 

mates  for  P,d,m  were  not  'good'  initial  estimates  for  the  Newton  proce¬ 
dure  they  are  listed  below  lor  completeness.  The  chi-square  values  are 
given  in  Table  V. 


Bus 

P-1.42526 

0*. 66959 

iw»2. 53973 

CDRT 

P—  .612943 

<*.84179 

mp3. 85879 

OEK 

P—1.77332 

0*. 94971 

m- 18. 76879 

HUO 

P-1.38298 

<*.86383 

m- 1.20253 

I  MFK 

P— .85582 

<*.98331 

mpl.81454 

MPDG 

P—1.94282 

<*.98172 

m»2. 17709 

MPDS 

P—2.25281 

<*1.82966 

m- 1.64836 

MTU 

P— 1.26428 

d* 1.18823 

m«4. 74657 

Processors 

P—1.29783 

<*1.18881 

np5. 73843 

SCU 

P—1.83275 

d«. 87388 

mp 1.84555 

SLU 

P— .91822 

<*1.15546 

m=2. 54579 

Table  V 


-Souare  Tail  Probabilities  -for  the 


Divergent  Three  Parameter  M 


Value 


[:<1 


For  the  items  which  did  converge,  using  the  Newton  procedure,  there 

s  A 

is  a  clear  trend  For  the  values  of  m  to  be  close  to  8  and  for  P  to  be 

a  a 

around  3.  But,  for  those  items  which  diverged,  m  is  larger  than  1  and  P 

is  very  small.  (Remember  P  is  actually  an  estimate  of  the  dose  which  in 

our  case  is  the  Log  WT8F)  . 

Comparing  Table  IV  with  Tables  I  and  II  it  can  be  seen  that  only 
one  item,  the  DSMU,  had  a  better  fit  under  the  three  parameter  model 
than  either  logit  or  probit  provided.  Two  other  items,  the  MFK  and  MHP, 
had  better  fits  than  with  logit.  All  fits  were  worse  than  the  probit 
except  the  DSMU. 

Looking  at  Table  V  shows  that  all  of  the  items  had  worse  fits  than 
that  provided  by  the  logit,  and  most  had  worse  fits  than  the  probit. 

But,  these  items  were  the  divergent  ones  and  their  fits  are  questionable 
anyway. 


Conclusions  and  Remarks 

The  technique  offered  by  Prentice  is  interesting  to  say  the  least. 
It  contains  within  its  -family  both  the  logistic  and  normal  distributions 
which  are  so  o-ften  used.  Prentice  offers  a  standard  score  test  to  exam¬ 
ine  the  fit  of  the  logistic  and  normal  models  which  is  more  sensitive 
than  the  usual  chi-square  test  of  fit.  This  test  is  based  on  the  asymp¬ 
totic  distribution  of  £'*(d_|/dm,djl/dn>  evaluated  at  (m,n) . 

However  sensitive  the  above  test  may  be  the  computations  for  it 
depend  on  convergence  of  the  Newton  procedure  so  that  the  information 
matrix  may  be  obtained.  Yo’<  must  also  reparameterize  as  m  and  n  ap¬ 
proach  infinity  for  the  normal  model.  But,  I  must  adnit,  this  is  not 
really  a  drawback  if  you  have  good  estimates  of  P  and  d. 

The  application  of  the  three  parameter  model  did  not  give  better 
results  than  those  obtained  by  the  more  familiar  probit  and  logit  meth¬ 
ods,  except  for  a  few  items  noted  previously.  Even  for  these  few  items 
the  computational  techniques  were  extremely  more  complex  and  time  con¬ 
suming.  However,  this  does  not  imply  that  the  technique  must  always  be 
so  diff icul t . 


I  feel  that  a  future  effort,  dedicated  to  computerizing  the  entire 

family,  would  be  useful.  There  were  many  numerical  techniques  (Refs  23, 

25,27)  which  I  did  not  investigate,  due  to  time  considerations,  which 

may  have  been  useful  in  overcoming  the  computational  difficulties  which 
I  encountered. 


Copen haver  and  Mielke  (Ref  11)  give  another  family  of  distribu¬ 
tions,  analogous  to  that  of  Prentice  (Ref  30).  This  family  is  desig¬ 
nated  as  the  omega  distribution.  The  cdf  of  this  distribution  is 

F(x(q)>  «  q  (5.1 


its  pdf  is 


f  (x(q)  )  »  1  -  |2q  -  1 1 


v+1 


and 


■A 


.-1 


x<q)  a  /  [f (x(z) > 3  4  dz 


(5.2 


(5.3 


where  8<q<l  and  v>-l.  As  was  the  case  for  the  family  given  by  Pren¬ 
tice,  this  family  has  embedded  in  it  special  cases.  "In  particular, 
this  distribution  is  a  double  exponential  when  va0,  a  logistic  distribu 
tion  when  v»l,  and  a  uniform  in  the  limiting  case  as  v-h»*  (Ref  11:177). 
That  is,  for  v»l  the  pdf  is  f (x(q))*4q( 1-q)  ,  and 

x(q)  » J [4z(l-z)J  1  dz  *  ( 1/4) log(q/l-q) 

•'va 

or  the  logit  of  q.  This  yields  the  logistics  density  function: 

f(x)  ■  4  exp(4x)  Cl  ♦  exp(4x)l  ^ 


For  v«8  the  double  exponential  density  function  is  f (x)=exp(-2 |x  |) , 
and  for  v-h*  the  density  is  f(x)»l  for  -l/2Cx<l/2. 


Note  that  the  symmetric  distributions  of  Eq(4.2)  (e.g  when  m=n) 
given  by  Prentice  are  approximately  the  '...subset  of  omega  distribu¬ 
tions  where  0<v<2’  (Ref  11:178).  But,  as  can  been  seen  by  Eq(5.2) ,  the 
omega  distribution  includes  no  asymmetric  distributions. 

The  likelihood  function  for  quantit  analysis  is  the  same  as  that  in 
Eq(4.3>,  but  now 

P.  =  /  f(t)  dt  =  F(or+£x.) 

1  •'.«  1 

Thus,  the  tolerance  distribution  is  given  by 

fta+fcr)  *  1-|2P.-1|V+1 

and 

a *fht.  *  h  (P.)  »/  l/(l-|2z-l|v+1>  dz  (5.4) 

1  v  i  •'i/2 

where  hv(P.)  i*  termed  the  'quantit"  of  Pi . 

Computations  for  the  Model 

"The  computational  procedures  for  obtaining  ML  estimates  of  param- 

A  A  yv 

eters  a,*  and  v  (a,*,v>  consists  of  an  efficient  search  routine  for  de- 

A 

termining  v"  (Ref  11:178).  The  calculations  are  identical  to  that  for 
Prentice  (see  Section  IV),  however;  the  procedure  is  more  tracktable 
since  the  numerical  calculations  are  simpler,  due  to  "nice'  functions. 
The  procedure  starts  at  v°"l  and  continues  until  0  is  of  the  de- 

A 

sired  accuracy  or  v  exceeds  20.  "If  v  is  20  or  larger,  there  is  very 
little  difference  between  the  omega  distribution  and  limiting  uniform 


distribution. . . 


(Ref  11:178).  For  the  initial  value  of  v<*»l  use  as 


initial  estimates  of  a  and  P  the  least  squares  solution  of  the  equation 


j 

r. 


h^CP)  ■  o+Px,  where  hv<P>  is  the  observed  quantit  corresponding  to  the 
proportion  of  observed  responses  s,  at  dose  x,  to  the  number  of  trials 
at  dose  x. 

The  derivatives  of  the  log-likelihood  and  the  elements  of  the  in¬ 
formation  matrix  are  calculated  exactly  as  in  equations  4.5  and  4.6  re¬ 
spectively.  However,  for  the  omega  model  with  P.»F(a+Px . )  the  deriv¬ 
atives  with  respect  to  the  parameters  a,P  are 

dP/da  *  f<a+Px) 


I 

!: 

i  e 

s 

C 

c 

I 


and 


dP/dP  -  xf  (ct+Px) 


So,  if  you  have  the  jth  iterative  estimates  of  a. and  P.and  (j-l)st  it- 

J  J 

erative  estimate  of  then  the  Newton-Raphson  procedure  yields  the  jth 

estimate  of  P. .  The  <m+l)st  iterative  solution  of  P.  .  is  given  by 
i  i ,  J 

p<m*l)  =  p<m>  +  ( 1-|2P^~  1|V+1  >  [<x.+  P . x .  —  h  <P<m)>] 
xj  U  xj  j  j  i  v  u 

and  since  we  can  write  Eq<5.4>  as 


i 


K 


and 


hence 


Q<P.  .)  ■  a.+P.x  -  h  <P.  .)  »  8 
xj  j  J  i  v  xj 


G'(P.  .)  *  -(1  -  i2P.  .-  l|v+i>“1 

1  J  1  J 


pfm+1)  =  p. .-0<p?">)/0'<p??)) 
i  j  xj  u  xj 
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Since  the  functions,  and  thus  the  above  equations,  are  'well-behaved* 
the  2  by  2  information  matrix  never  had  the  problems  that  I  encountered 
with  Prentice's  model. 

The  only  computational  difficulty  arises  from  the  need  to  calculate 
the  quantit  h^(P) .  Copen haver  and  Mielke  initially  used  the  infinite 

series  given  by 


where 


h  <P>  * 
v 


£Y_i2e±iii<vM>M 

2  zl  i<v+ 1)  + 1 

i«e 


6 


1  P> 1/2 
e  p«i/2 
-1  P<l/2 


"Whenever  v>l  and  |2P-1|<0.9,  the  ...  series  converges  rapidly.  In 
fact...  the  maximum  error  will  be  10-*  when  the  first  30  terms  are  sum¬ 
med*  (Ref  11:186).  However,  when  v<l  or  |2p-l|  is  close  to  one  (as  in 
our  case)  you  sum  for  an  'appropriate'  number  of  terms  and  then  add  a 
remainder.  The  remainder  is  in  the  form  of  a  continued  fraction.  "The 
continued  fraction  converges  slowly  when  either  v  is  near  -1  or  |2p-l| 
is  close  to  one"  (Ref  28:222).  Magnus  et.  al.  (Ref  28)  develop  a  closed 
form  expression  for  the  sum  of  an  infinite  series.  Copenhaver  and  Miel¬ 
ke  adopted  this  latter  technique  in  their  computerization  (Appendix  E 
contains  the  source  code  as  written  by  Copenhaver  and  Mielke). 


Results  of  the  Quantit  Analysis 


The  values  obtained  for  a,J*  and  v  are  given  below  for  the  17 


equipment  items 

f<cc+4x)  is  the 

,  where  or+Fx.  ■  h  <P> 

I  V 

tolerance  distribution 

as  given  by  Eq<5.4) 

e 

,  and 

ART 

v-19.08 

a— .37446 

4— .02797 

Bus 

v*28 

OF-. 33543 

4—. 03568 

CDRT 

v— .90 

of9. 05674 

4—5.91737 

DEK 

v- .90 

cf8. 49762 

4— 5.23146 

DSMU 

v«29 

of-. 37875 

A— .02647 

HUD 

v«20 

of-. 35067 

4—. 02975 

IMFK 

v— .90 

of  1 . 69570 

4—5.96704 

INS 

v=28 

of-. 37296 

4—.  02658 

MFK 

v»20 

OF-. 40607 

4—. 02475 

MiP 

v»20 

OF-. 39221 

4— .02368 

MPDG 

v=20 

OF-. 442 13 

4— .02624 

MPDS 

v=*20 

a— .44539 

4— .02371 

MTU 

v=-.90 

0F8.45196 

4— 4.19618 

Processors  v— .98 

CF8.44919 

4—4.48845 

SCU 

v»20 

of-. 355 13 

4— .82874 

SLU 

v— .90 

cf6. 63052 

4—4.30218 

SMRT 

v*»9.50 

of-. 38000 

4— .05279 

The  value  of  v  for  all  of  the  above,  with  one  exception  the  SMRT, 
is  either  -.9  or  28.  These  results  are  somewhat  surprising,  since  fcr 
those  items  with  v— .9  the  tails  are  very  heavy.  The  clear  implica¬ 
tion  is  this;  an  item  with  a  parameter  value  of  -.9  for  v  has  a  very 
narrow  band  of  critical  MTBFs.  For  MTBFs  below  this  narrow  band  the 


item  will  always  cause  an  abort  and  'for  MTBFs  above  it  the  item  will 
never  cause  an  abort.  Considering  that  the  proportion  of  aborts  is  very 
low  in  the  sample  data,  and  somewhat  linear  (flat),  this  is  surprising 
indeed. 

Now  looking  at  Table  VI  and  comparing  these  chi-square  values  with 
those  in  Tables  I  and  II  for  probit  and  logit  respectively,  it  is  clear 
that  the  quantit  method  gave  consistently  better  fits  than  the  probit. 

Hi th  only  two  exceptions,  the  MPDG  and  the  MTU,  quantit  also  gave  better 
fits  than  did  logit. 


Table  VI 


Chi -Sou are  Tail  Probabilities 
for  Quantit 


Itsm 

Value 

IMFK 

.9915 

CDRT 

.9846 

Bus 

.9778 

MFK 

.9162 

SMRT 

.8788 

Processors 

.3781 

MiP 

.8271 

SLU 

.5994 

ART 

.5738 

MPDG 

.4920 

DEK 

.4499 

MPDS 

.3994 

HUD 

.2111 

DSMU 

.1895 

INS 

.1680 

SCU 

.0815 

MTU 

.0255 

Conclusions  and  Remarks 

This  method  of  quanta!  analysis  appears  to  have  been  better  than 
all  of  the  preceding  methods.  It  gave  consistently  better  fits  than 


probit,  one-hit,  and  the  generalization  of  Prentice  for  all  1?  equipment 
items.  It  also  had  better  fits  for  15  of  the  1?  items  than  did  logit. 

I  found  that  quantit  analysis  was  quite  easy  to  implement  and  un¬ 
derstand.  The  computer  code,  given  to  me  by  Hr.  Copenhaver,  would  be 
'easy'  to  modify  to  attain  any  accuracy  of  the  parameters  needed.  It 
should  be  noted  that  the  code  contains  a  subroutine  to  perform  probit 
analysis  (logit  is  a  by-product  of  the  program;  v*l  is  the  logit  and  is 
always  output) . 

I  would  recommend  that  a  future  effort  consider  modifing  the  code 
so  that  instead  of  doing  fixed  iterations  on  v,  it  performs  a  simulta¬ 
neous  search  for  all  three  parameters  a,P  and  v.  However,  this  would 
then  have  the  same  sensitivity  problems  as  Prentice's  model.  I  would 
further  recommend  that  Prentice's  model  be  incorporated  with  this  code, 
thus;  giving  a  more  comprehensive  analysis  tool. 


VI  Two  Families  of  Trans-formations 


Previously  the  logistic  model  has  been  a  special  case  of  all  the 
other  models.  As  a  matter  of  -fact,  the  determination  of  how  well  these 
other  models  performed  has  been  a  comparison  of  them  against  the  logis¬ 
tic,  via  chi-square  tail  probabilities.  However,  the  logistic  model  is 
only  a  tentative  model  and  we  need  to  look  at  how  adequate  it  is.  "If 
we  can  find  a  procedure  which  detects  inadequacy,  and  also  indicates  the 
kind  of  desirable  modification  to  the  mode),  this  is  potentially  useful" 
(Ref  1:357). 

Mr.  Aranda-Ordaz  (Ref  1)  gives  us  two  families  of  models  which  help 
achieve  the  stated  objective.  These  families  each  contain  the  logistic 
and  alternatives  as  special  cases.  These  models  each  have  an  associated 
transformation  which  model  symmetric  and  asymmetric  departures  respec¬ 
tively  from  the  logistic.  Symmetric  transformations  are  such  that  they 
lead  to  "...essentially  the  same  answers  if  successes  and  failures  are 
interchanged"  (Ref  1:357). 

The  Symmetric  Family  and  Associated  Test 

The  symmetric  family  of  alternatives  to  the  logistic  is  given  by 

t  (e>  =  (2/x)/ex-(  1-0) x\  (<s.i> 

\0X*(1-0)X/ 

where  0  is  the  probability  of  success  and  x  is  the  transformation  para¬ 
meter.  If  we  let  T^<0)=A+B  Ln  x  and  then  solve  Eq(6.1)  for  0  as  a 


function  of  t  we  then  obtain 


9<t) 


<  l+Xf/2) 


1/x 


< l+XT/2) 1/X+( 1-XT/2) 1/X 


V  1 


(xt/2<-i) 


|Xt/2 I  <1  (4.2) 


(Xf/2*1> 


Note  that  Eq(4.1)  '...reduces  to  the  logistic  transformation  in  the 
limit  when  X=0  and  to  the  linear  transformation  when  X=»l"  (Ref  1:358). 
Another  feature  of  Eq(4.1)  is 


T  <$)  *  -T  (1-0)  ;  T  (6)  =  T  ($> 

A  o  A  — A 

or  the  treatment  of  successes  and  failures  is  symmetric. 

As  is  the  case  for  all  quantal  methods  the  underlying  process  is 
binomial  and  the  likelihood  function  is  the  same  as  that  in  Eq(4.4) . 
The  systematic  part  of  the  model  is  given  by 


If  we  fit  by  maximum  likelihood  a  linear  expression 
for  T  for  a  range  of  values  of  X,  we  may  consider 
the  maximized  log  likelihood  as  a  function  of  X  and 
hence  derive  not  only  the  maximum  likelihood 
estimate  x,  but  also  determine  which  values  of  X 
provide  an  acceptable  fit.  (Ref  1:358) 


Since  it  was  presumptuous  to  assume  the  logistic  as  the  model  it  is 
just  as  presumptuous  to  assume  the  symmetric  will  yield,  via  MLE,  any¬ 
thing  but  the  logistic  model  back.  Hence,  before  proceding  with  maximum 
likelihood  estimation  of  x,  we  need  some  test  to  determine  if  there  are 
indeed  symmetric  departures  from  the  logistic  model.  The  hypothesis  we 
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want  to  test  is  H:  x=8. 


If  we  consider  the  parameter  vector  £  as  a  nuisance  parameter,  and 
replace  it  by  its  maximum  likelihood  estimate  under  X=0,  then  the  test 
statistic  is  the  efficient  score  U(X)*dl/dX.  U(x>  vanishes  at  x=0,  but 
after  some  reparameterizing,  the  score  takes  the  following  limiting 
form: 

U<0)  =  £(r.-n.0.)t»/12  (4.3) 

isl 

A 

where  0.=l/( l+exp(-T.))  ,  T.  is  the  logit  equation,  r.is  the  number  of 
responses,  and  ni  is  the  number  of  trials. 

Since  Eq(6.3)  is  distributed  asymptotically  normal  (Ref  14:363)  the 
test  may  be  carried  out  with  the  standardized  form  of  Eq(6.3) .  The  re¬ 
jection  of  H:  X=0  is  for  large  values  of  the  test  statistic.  The 
variance  of  Eq(6.3)  is  given  by 


I  (0)  = 

n  xx  pP  0x 


(6.4) 


where 


and  1^  are  the  partition  submatrices  of  I, 


Fisher's 


information  matrix  (see  Section  IV,  in  particular  Eq(4.6>  for  more,  or 
DeGroot  (Ref  14)  section  7.8  ).  The  individual  components  are 


isl 

m 


( s—  1 , . .  •  ,p) 


where  d.=  0.<l-0.),  and  p  is  the  dimension  of  p.  In  our  case  this  di- 
mension  is  2  and  the  Xj  vectors  are  the  unit  vectors.  Note,  "All  values 
required  to  perform  the  test  may  be  computed  from  the  output  of  a  logis¬ 
tic  fit"  (Ref  1:340) . 


The  Asymmetric  Family  and  Associated  Test 


lihood,  again  just  as  -tor  the  symmetric.  We  also  develop  a  test  statis¬ 
tic  analogous  to  that  of  the  symmetric.  The  hypothesis  is  now,  however, 
Hs  \>1.  The  test  statistic  is  given  by 


m 


U< 


1)  »\  <r.-n.0.)  <0.-Hog<  1-©.)) 
)  i  i  x  i  i 

L-i  0. 
i  =  l  1 


(6.7) 


where  exptTJ/t  l+exp<T.)  ]  . 

The  variance  is  given  by  Eq(6.4>,  but  now  the  individual  components 

are 


I  a  £n.<0.  +  log<l-0.>>2/exp<T.> 

XX  11  1  1 

m 

1  =  £<0.*log<l-0.>)n.x  <1-0. >  <s=l,...,p> 

s  «=1  1  i  x  s  l 

lft  p  9  EVsW^V  <r’s^1 . p> 

s  P  ill 

again  p  is  the  dimension  of  £.  This  test  also  requires  only  values  from 
the  logistic  fit.  The  rejection  of  the  hypothesis  is  for  large  negative 
values  of  the  standardized  statistic. 

Computational  Procedures 

The  calculations  for  both  the  test  statistic  UOJ  and  the  two  fami¬ 
lies  are  straightforward.  The  calculations  for  U(x>  where  explicitly 
given  in  the  previous  paragraphs  and  no  further  explanation  will  be  giv¬ 
en  here,  except  for  one  small  discussion. 

It  is  common  practice  to  use  weighted  feast  squares  procedures  when 
performing  logit.  While  Mr.  Arand-Ordaz  does  not  explicitly  state  that 
he  uses  weighted  least  squares  I  will  assume  that  he  did.  The  reason 
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for  weighting  in  a  least  squares  procedure  is  that  the  assumption  of  ho- 
moscedastici ty  (constant  variance)  may  not  be  justified  (Cox  (Ref  12) 
section  2.2  and  Thei I  (Ref  34)  section  6.2).  The  assumption  of  homosce- 
dasticity  is  one  of  the  key  assumptions  underlying  all  regression  analy¬ 
sis,  and  must  not  be  ignored.  However,  for  all  quantal  analysis,  since 
the  underlying  process  is  binomial,  weights  are  easily  obtainable.  Since 
the  variance  of  a  binomial  distribution  is  npq,  it  makes  sense  to  use 
this  as  the  weighting  coefficient,  which  is  what  I  have  done.  This  is 
also  what  Dr.  8arr  (Ref  3)  has  done  and  what  is  suggested  by  Finney  (Ref 
15)  . 

The  model  calculations,  as  1  said  before,  are  straightforward  once 
you  have  determined  if  there  are  symmetric  or  asymmetric  departures  from 
the  logistic.  First  for  a  fixed  X  calculate  the  value  for  Eq  (6.1)  or 
Eq  (6.5),  for  each  level  of  the  stimulus  (MTBF) .  Next  using  these  val¬ 
ues  perform  a  least  squares  regression  to  obtain  ^  and  the  regres¬ 
sion  coefficients.  Then,  using  this  regression  equation,  calculate  6(T) 
in  Eq(6.2>  or  Eq(6.6) .  Finally,  using  the  values  of  9<T)  just  calculat¬ 
ed,  determine  the  value  of  the  log  likelihood  function.  Since  you  are 
trying  to  maximize  the  likelihood  function  the  value  of  x  that  yields 
the  largest  value  of  the  likelihood  is  the  MLE  estimate  x. 

Note,  one  would  normally  only  apply  one  of  the  two  families  if  the 
test  statistics  indicated  a  rejection  of  the  null  hypothesis  that  the 
distribution  is  logistic.  However,  for  comparison  purposes  I  calculat¬ 
ed  the  MLE  for  all  17  equipment  items  for  both  families. 

Results  for  the  Symmetric  Family 


Given  in  Table  VII  are  the  results  of  the  normalized  test  statistic 


computed  -from  Eqs  (6.3)  and  (6.4),  and  listed  in  ascending  order 
listed  are  the  MLE  values  for  X. 


Table  VII 


Item 

Statistic 

X 

IMFK 

.5885 

0. 

CDRT 

2.0299 

0. 

DEK 

2.4989 

0. 

SMRT 

2.8950 

0.13 

Processors 

2.9858 

0. 

MFK 

4.1348 

0.32 

Bus 

5.0072 

0. 

ART 

5.1847 

8.18 

SLU 

4.0953 

0. 

INS 

6. 1689 

0.22 

MMP 

7.2325 

0.33 

DSMU 

10.3723 

0.32 

MPDS 

13.4003 

0.24 

HUD 

14.4408 

0.62 

MPDG 

15.6428 

0.33 

SCU 

29.4729 

0.88 

MTU 

35.8081 

0. 

Comparing  the  order  of  Table  VII  and  Table  II  (logit)  it  is  inter¬ 
esting  to  note  that  the  order  is  somewhat  the  same,  particularly  for  the 
tops  and  bottoms  of  the  tables.  This  is  what  one  would  expect  since,  as 
the  logistic  fit  becomes  worse,  the  test  statistic  should  indicate  de¬ 
partures  from  the  logistic. 

Remember  the  null  hypothesis  is  H:  x=8,  or  the  distribution  is  not 
different  from  the  logistic.  Therefore,  the  test  indicates  whether  or 
not  there  are  symmetric  departures  from  the  logistic.  The  values  of  x, 
in  Table  VII,  are  fairly  consistent  with  the  expected  results  as  indi¬ 


cated  by  the  test  statistic. 

I  am  sure  you  have  already  noticed  that  the  values  of  the  test  sta- 


tistic  are  larger  than  one  would  expect  -for  a  standard  normal  distribu¬ 
tion.  This  could  be  due  to  the  underlying  distribution  not  being  logis¬ 
tic  or  even  shaped  liked  a  one,  but  something  completely  unrelated  to 
the  logistic.  Indeed,  the  underlying  distribution  may  not  be  represen¬ 
ted  by  the  -family  in  Eq<6.1)  at  all.  I-f  this  is  the  case  the  test  for 
departures  is  not  valid.  This  may  also  explain  the  zero  values  for  x 
for  those  items  with  very  large  test  values. 

Listed  below  are  the  equations  for  T,  attained  at  the  respective  x 
value,  for  all  17  items. 


© 


IMFK 

T= 

CDRT 

T= 

DEK 

T= 

SMRT 

T= 

Processors 

T* 

MFK 

T=- 

Bus 

T= 

ART 

T=  ■ 

SLU 

T= 

INS 

T«  • 

M1P 

t=-: 

DSMU 

t»- 

MPDS 

T— 

HUD 

T»- 

MPDG 

T=- 

SCU 

T=- 

MTU 

T«  ■ 

.88418  -  1.08979  In  x 
.49617  -  1.16095  In  x 
.34701  -  1.01782  In  x 
.86791  -  .91811  In  x 
.50983  -  .89628  In  x 
.98938  -  .44583  In  x 
.13606  -  1.00295  In  x 
.39222  -  .78307  In  x 
.00023  -  .84037  In  x 
.57526  -  .69058  In  x 
.02101  -  .37681  In  x 
.48325  -  .47681  In  x 
.33997  -  .40853  In  x 
.18344  -  .13002  In  x 
.02101  -  .37681  In  x 
.12457  -  .01772  In  x 
.74551  -  .66784  In  x 


Listed  in  Table  VIII  are  the  chi-square  values  for  the  fit  of  the 


17  equipment  items  using  the  symmetric  family.  The  values  are  listed  in 
descending  order.  Comparing  the  values  of  Table  VIII  with  those  of  Ta¬ 


ble  II  the  order  is  essentially  the  same.  However,  there  is  an  increase 
in  the  value  of  the  tail  probabilities  in  Table  VIII  for  those  items 
which  did  not  have  an  MLE  \=0  (remember  if  X»0  these  items  are  fit  by 
the  logistic,  hence,  the  same  as  Table  II>.  There  are  three  exceptions 
however.  These  exceptions  are  the  SMRT,  ART,  and  INS,  but  their  tail 
values  (Table  VIII)  are  almost  identical  with  those  for  logit  in  Table 


Table  VIII 


I  tern 

IMFK 

CDRT 

Bus 

MFK 

MPDG 

l“MP 

SMRT 

Processors 

MPOS 

ART 

SLU 

DEK 

SCU 

HUD 

DSMU 

INS 

MTU 


Value 

.9754 

.9728 

.9641 

.9389 

.9075 

.8261 

.8247 

.8241 

.6873 

.4592 

.4527 

.3617 

.3107 

.2148 

.1866 

.0799 

.0487 


Table  IX  lists  the  results  of  the  normalized  test  statistic  com- 
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puted  from  Eqs<6.7)  and  (6.4),  and  are  in  ascending  order.  Also  listed 
are  the  values  of  the  MLE  x. 

Table  IX 

Values  for  the  Asymmetric 


1 


a 


Item 

Statistic 

X 

DEK 

-1.2765 

-138 

CDRT 

-  .7509 

-148 

Processors 

-  .5034 

-8 

IMFK 

-  .2972 

-267 

Bus 

-  .0467 

-781 

SMRT 

.1812 

20 

SLU 

.1970 

-147 

MFK 

.3151 

227 

MPDG 

.9758 

-1789 

HUD 

1.1261 

510 

MTU 

1.3603 

-531 

MMP 

1.4385 

325 

ART 

1.5993 

61 

INS 

3.1261 

59 

SCU 

3.4055 

3490 

DSMU 

3.7150 

223 

MPDS 

4.0440 

859 

z 

i 


2 

i 


The  theoretical  null  hypothesis  in  this  case  is  H:  x*l  ,or  the 
logistic.  Therefore,  for  large  values  of  the  test  statistic  you  will 
reject  the  null  hypothesis.  However,  Hr.  Aranda-Ordaz  indicates  that 
you  will  reject  only  for  large  negative  values.  This  implies  that  the 
working  hypothesis  is  really  Hs  Xil.  The  reason  for  this  is  clear; 
since  we  are  interested  in  asymmetric  departures  from  the  logistic  we 
are  really  interested  in  values  of  X  significantly  less  than  one.  While 
values  larger  than  one  are  also  acceptable  the  distribution  starts  to 
lose  its  asymmetry  and  becomes  more  symmetric  looking. 

If,  as  I  said  before,  I  were  to  go  by  the  test  statistic  value,  as 
an  indication  whether  or  not  to  proceed  with  the  asymmetric  model,  I 
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would  not  proceed.  I  would  be  hard  pressed  to  reject  the  null  hypoth¬ 
esis  H:x>l  for  any  of  the  items,  except  possibly  the  DEK.  But,  since  we 
are  interested  in  how  well  the  over  all  asymmetric  model  and  test 


perform  for  our  data  1  continued,  and  obtained  the  following  equations 
for  T.  Again  these  equations  are  attained  at  the  MLE  X  given  in  Table 
IX. 


DEK 

T*- 1.74791 

— 

.72424 

In 

CDRT 

T*  -.99838 

- 

.94095 

In 

Processors 

T*  .22093 

- 

.85761 

1  n 

IMFK 

T*-2. 48298 

- 

.77771 

1  n 

Bus 

f»-2. 24845 

- 

.75122 

1  n 

SMRT 

.93780 

- 

1.09026 

In 

SLU 

1 .7 1897 

- 

.62746 

In 

MFK 

T=*  3.01833 

- 

1.45962 

In 

MPDG 

T^*-5 . 4457 1 

- 

.45658 

In 

HUD 

T*  12. 85090 

- 

2.66024 

In 

MTU 

1^-3. 93047 

- 

.34990 

In 

MMP 

T»  3.97039 

- 

1 .39584 

1  n 

ART 

T«  1.46151 

- 

1 . 16937 

In 

INS 

T®  2.02577 

- 

1 . 17478 

In 

SCU 

T®54. 15579 

- 

7.87826 

In 

DSMU 

T-  5.40580 

- 

1.67762 

In 

MPDS 

T-  .34738 

— 

1 . 17886 

In 

x 

x 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


Table  X  contains  the  chi-square  tail  probabilities  for  the  fit  of 


the  data  using  the  asymmetric  'family.  As  with  all  previous  chi-square 
tables  these  also  are  listed  in  descending  order.  Note,  the  order  is 
almost  identical  to  that  of  Table  Mill  and  also  Table  II.  With  four  ex¬ 
ceptions  these  tail  probilities  are  better  than  those  in  Table  Mill, 
hence,  better  than  those  of  Table  II.  The  four  exceptions  are  the  MFK, 
SMRT,  SCU,  and  HUO. 

Comparing  the  tail  values  for  these  four  items  with  those  in  Table 
Mill  shows  that  while  the  symmetric  was  better  it  was  not  significantly 
better.  The  only  item  of  the  four  which  was  not  relatively  the  same  as 
the  logit  fit,  in  Table  II,  is  the  SCU.  For  the  SCU  both  the  symmetric 
and  asymmetric  models  gave  chi-square  values  more  than  four  times  larger 
than  that  of  logit,  a  marked  improvement. 


Table  X 


for  the  Asymmetric  Famil 


Processors 


Malue 


One  last  interesting  note  concerns  the  test.  If  I  had  performed 


the  asymmetric  test,  with  the  intention  o f  only  applying  the  model  if 
indicated  by  the  test,  the  only  item  I  would  have  tried  to  fit  would 
have  been  the  DEK.  And,  for  this  single  item,  I  would  have  found  an 
almost  two-fold  increase  in  the  chi-square  tail  probability.  But,  I 
would  have  missed  improved  fits  for  many  other  items. 

Conclusions  and  Remarks 

I  feel  the  numerical  calculations  which  I  performed  are  correct  and 
accurate,  for  both  the  test  statistic  values  and  the  models  themselves. 
However  I  would  be  less  than  honest  if  I  did  not  report  the  following 
discrepancy. 

Mr.  Arand-Ordaz  applied  the  asymmetric  test  to  one  set  of  sample 
data.  This  data  has  a  point  where  the  responses  are  the  entire  sample 
(e.g.  the  number  of  insects  killed  equaled  the  number  exposed  to  the 
chemical  stimulus)  (Ref  1:342  table  4).  How  he  treats  this  point  was 
not  indicated.  This  point  is  a  problem,  as  a  matter  of  fact  any  point 
where  all  or  none  of  the  subjects  respond  causes  problems.  There  are 
many  schools  of  thought  about  what  to  do  with  a  point  like  this,  but  the 
standard  seems  to  be  to  replace  p»s/n  by  p»(s+.5)/(n+l)  (Ref  11:178). 
However,  if  one  does  not  use  this  method  than  something  very  similar  to 
it  is  usually  suggested. 

With  the  above  in  mind  I  tried  several  different  methods,  but  could 
never  attain  his  fitted  values  as  stated  (Ref  1:342  table  4)  for  the  lo¬ 
gistic  model.  I  could  get  very  close  using  the  SPSS  regression  proce¬ 
dure  to  perform  the  logistic  fit,  but  never  close  enough  to  satisfy  my¬ 
self.  Normally  I  would  chalk  this  up  to  not  knowing  how  Mr.  Arand-Ordaz 
treats  the  problem  point.  However,  to  calculate  the  test  statistic  you 


need  the  logistic  -tit.  And,  I  could  not  duplicate  this  either,  at  least 
without  "cheating*  as  I  will  explain. 

Mr.  Aranda-Ordaz  (Ref  1:362)  reports  a  test  value  o-f  -2.76  for  the 
asymmetric  statistic.  Using  the  closest  logistic  fit  I  could  attain 
(compared  to  his)  I  could  not  get  this  value  (remember  this  value  is 
obtained  by  the  standardized  form  of  U(x) ,  which  has  a  mean  of  zero  and 
variance  as  given  in  equation  6. A.  However,  I  could  get  his  value  if  I 
divided  U(X)  not  by  the  standard  deviation,  but  by  the  variance.  Since 
the  equations  ( 6.7  and  6.4)  are  easy  to  calculate  and  verify,  and  since 
the  SPSS  regression  technique  is  valid,  I  am  convinced  that  the  value 
of  -2.76  is  wrong  as  reported  by  Mr.  Arand-Ordaz. 

These  two  sets  of  models  were  relatively  easy  to  apply  and  they  in¬ 
corporate  some  of  the  same  calculations  required  in  Sections  III  and  IV. 
Hence,  these  models  could  easily  be  incorporated  into  a  larger  computer 
package  which  included  the  families  given  by  Prentice  (  Ref:31)  and  Cop- 
enhaver  (Ref  11).  The  three  families  together  include  the  logistic  (lo¬ 
git),  normal  (probit),  and  numerous  symmetric  and  asymmetric  alterna¬ 
tives.  A  comprehensive  quantal  assay  computer  package  has  much  poten¬ 
tial  for  use  in  many  areas  other  than  bio-assay 


Appendix  A 


Description  of  the 
Avionics  Evaluation  Prooram 
AEP  Model 

The  Air-to-Ground  Mission  Analysis  (MAP)  submodel  of  The  AEP  evalu¬ 
ates  the  performance  of  a  flight  of  up  to  four  aircraft  on  a  mission 
which  may  involve  multiple  targets,  multiple  search  passes,  and  multiple 
attack  passes.  The  aircraft  proceed  along  an  externally  generated  nomi¬ 
nal  trajectory  through  the  mission  phases  of  takeoff,  navigation  to  the 
search  area,  search,  attack,  and  return  to  base.  Monte  Carlo  techniques 
are  applied  to  mean-time-between-fai lure  (MTBF)  data  for  the  defined  a- 
vionics  throughout  the  mission  to  determine  which  subsystem  modes  are 
functioning,  restoring  to  back-up  modes,  and  mission  aborts  as  required. 
Target  location  uncertainties  and  navigation  system  performance  parame¬ 
ters  are  combined  to  define  the  actual  flight  path  relative  to  the  true 
target  location.  The  sensor  ground  swath  for  the  defined  search  pattern 
is  then  compared  to  the  true  target  location  to  determine  if  the  target 
passes  through  the  sensor  ground  area  coverage.  Probabilities  of  de¬ 
tection,  target  kill,  and  aircraft  survival  are  sampled  to  determine 
which  mission  phases  are  successfully  completed.  The  model  utilizes  the 
best  mode  still  available  for  each  function  at  the  time  it  is  to  be  per¬ 
formed. 

The  MAP  includes  a  detailed  model  of  the  ground  turnaround  process. 
Preflight,  thruf light,  postflight  maintenance,  ordinance  arming/loading, 
refueling,  scheduling,  and  launch  are  modeled  in  terms  of  time  require- 


merits  and  event  uncertainties. 

Prior  to  launch  of  the  first  sortie  of  the  day,  the  model  sequenced 
each  aircraft  through  preflight  maintenance  and  ordinance  loading.  Dur¬ 
ing  the  maintenance  interval  equipment  items  are  repaired  and  *time-to- 
repair"  is  recorded.  Prior  to  launch  of  subsequent  sorties,  the  model 
sequences  each  aircraft  through  ordinance  de-arming,  thru-flight  mainte¬ 
nance,  refueling  and  ordinance  loading.  Repair  and  loading  times  are 
recorded.  Upon  completion  of  the  last  sortie  of  the  day,  the  model  se¬ 
quences  each  aircraft  through  ordinance  de-arming,  post-flight  mainte¬ 
nance,  and  refueling  for  the  next  day. 

The  launch  subfunction  represents  the  interval  between  engine  start 
to  takeoff.  At  this  time,  an  additional  equipment  check  is  made  to  de¬ 
termine  additional  failures.  To  determine  these  launches  the  sortie 
scheduling  algorithms  utilizes  user  supplied  data  to  manage  the  starting 
time  for  individual  sorties  on  sequential  days. 

The  preflight,  thruf light,  and  postflight  maintenance  times  are 
based  on  mean  duration  time  input  data.  Ordinance  loading,  arming  and 
de-arming  times  are  determined  in  a  similar  manner.  Refueling  times  are 
based  upon  an  input  refueling  rate  (Ibs/min)  and  aircraft  fuel  storage 
capacity.  The  model  calculates  fuel  utilization  for  each  sortie  to  de¬ 
termine  additional  fuel  requirements. 

In  addition  to  ground  preparation  and  ground  maintenance,  the  user 
must  also  define  the  in-flight  mission  functions  along  with  their  var¬ 
ious  parameters  (e.g.  nav  accuracy,  drift  rates,  ect.)  and  their  associ¬ 
ated  suites  of  hardware  and  the  hardware  reliability  and  maintainability 
parameters.  The  in-flight  functions  are  navigation,  navigation  update, 
communication,  target  acquisition,  weapon  delivery,  general  flight,  tar- 


get,  end  survivability.  Note  that  each  of  these  functions  have  numerous 
subfunctions  (e.g.  navigation  has  radio-aided  nav,  and  self-contained 
nav  subfunctions). 

Once  airborne,  an  aircraft  must  have  one  of  two  navigation  func¬ 
tions  working  and  the  communications  function  working  or  the  aircraft 
will  abort.  The  other  functions  which  will  cause  an  abort  are  the  tar¬ 
get  acquisition,  weapon  delivery,  and  some  of  the  general  in-flight 
functions.  This  latter  function  is  one  that  can  be  used  to  determine 
additional  abort  conditions. 


Abort  Logic 


For  the  single  Aircraft  mission  simulation  the  concept  of  mode  re¬ 
gression  and  the  abort  process  is  straight  forward.  For  each  mode  with¬ 
in  each  function,  two  things  must  be  defined;  operating  or  performance 
characteristics,  and  a  suite  of  hardware  items  (which  also  has  a  set  of 
parameters  reliability  and  maintainability)  needed  to  perform  that  oper¬ 
ation  (mode)  within  the  function. 

Consider  one  of  the  functions  that  will  cause  an  abort  if  all  modes 
fail.  This  function  has  say  18  modes  and  within  each  mode  a  suite  of 
hardware  items  is  defined.  If  one  of  the  hardware  items  fails  in  mode  1 
then  the  aircraft  regresses  (moves)  to  mode  2;  you  can  think  of  this  as 
either  a  backup  mode  (e.g.  redundant  aircraft  equipment  items  may  have 
been  defined  in  this  mode),  or  a  degraded  mode  with  degraded  perform¬ 
ance  characteristics.  Now  if  a  crucial  hardware  item  fails,  that  is  one 
that  is  needed  by  all  modes,  or  if  enough  hardware  items  fail  such  that 
you  have  regressed  through  all  modes,  then  an  abort  will  occur.  Mode 
regression  starts  at  the  best  possible  mode  and  regresses  to  the  less 


desirable  modes. 
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Appendix 

C 

SamDle  Data 

o| 

The  -following  is  a  list  of  the  17 

equipment  items  considered  for 

analysis,  and  the  raw  data  generated  by 

the  AEP. 

i  u  * 

“S 

ART  -  Avionics 

Remote  Terminal 

a 

MTBF 

Launches  AttemDted 

Total  Aborts 

3888 

5378 

2 

2000 

5354 

6 

! 

1580 

5378 

3 

_ 

1000 

5437 

10 

. 

759 

5252 

12 

H 

694 

5340 

9 

r 

508 

5329 

17 

K< 

450 

5402 

13 

rl 

350 

5400 

24 

•*J 

250 

5316 

44 

i 

100 

5434 

78 

>' 

> 

w 

s: 

BUS  -  The  aircraft  data  bus 

> 

MTBF 

Launches  Attemoted 

Total  Aborts 

y 

© 

20800 

5443 

0 

1 

3000 

5380 

8 

'.- 

1500 

5309 

3 

V 

1800 

5386 

3 

\r 

• 

900 

5417 

6 

800 

5267 

6 

•  ' 

i 

700 

5277 

6 

■ 

600 

5275 

7 

>■ 

580 

5236 

10 

,V 

CDRT  -  Control 

and  Disci  ays  Remote  Terminal 

MTBF 

Launches  Attemoted 

Total  Aborts 

• 

4000 

5413 

1 

J* 

3800 

5396 

1 

i: 

2000 

5374 

2 

1500 

5390 

2 

1008 

5299 

1 

■ 

750 

5385 

3 

694 

5443 

5 

.• 

580 

5372 

6 

*•  * 

358 

5365 

11 

275 

5351 

13 

225 

5403 

13 

;a 

150 

5290 

23 

»*, 

• 

100 

5388 

46 

?■ 


Total  Aborts 


INS  -  Inertia)  Navigation  System 


O 


KTBF 

Launches  Attempted 

Total  Aborts 

2008 

5349 

8 

1500 

5444 

5 

1000 

5358 

13 

750 

5519 

8 

500 

5344 

21 

250 

5427 

47 

100 

5484 

115 

74 

5443 

110 

MFK  -  Multi 

-Function  Kfytjoar^J 

MTBF 

Launches  Attempted 

Total  Aborts 

3000 

5333 

1 

2080 

5472 

1 

750 

5240 

8 

540 

5329 

7 

350 

5301 

15 

100 

5204 

41 

75 

5315 

54 

MMP  -  Master  Mode  Panel 

MTBF 

Launches  Attempted 

Total  Aborts 

8808 

5381 

1 

4800 

5284 

2 

3000 

5358 

5 

2000 

5277 

2 

1508 

5372 

8 

1000 

5328 

14 

750 

5344 

15 

580 

5288 

22 

250 

5381 

34 

HPDG  -  Multi-Purpose  Display  Generator 


mue 

Launches  Attenuated 

Total  Aborts 

5800 

5309 

0 

3808 

5381 

0 

2008 

5341 

0 

1508 

5429 

0 

1120 

5443 

2 

900 

5248 

0 

750 

5298 

2 

540 

5274 

0 

450 

5341 

3 

350 

5299 

1 

100 

5311 

8 

44 


KTBF 

Launches  Attemoted 

Total  Aborts 

5800 

5356 

0 

3000 

5270 

0 

2000 

5412 

1 

1500 

5443 

0 

1000 

5386 

0 

750 

5385 

1 

500 

5416 

2 

350 

5264 

6 

250 

5333 

8 

100 

5329 

14 

75 

5226 

9 

J  -  Mul tiDlex 

Terminal  Unit 

HTBF 

Launches  Attenuated 

Total  Aborts 

17681 

5443 

7 

16000 

5389 

0 

8000 

5273 

3 

4000 

5383 

7 

2000 

5264 

10 

1500 

5237 

22 

1000 

5303 

27 

isors  -  Aircraft  data  processors 


HUE 

Launches  Attemoted 

Total  Aborts 

5000 

5240 

1 

4500 

5363 

6 

4000 

5272 

6 

3000 

5443 

10 

2000 

5365 

10 

1500 

5379 

11 

1250 

5433 

16 

1000 

5309 

18 

900 

5396 

16 

750 

5283 

18 

500 

5387 

35 

350 

5307 

47 

100 

5322 

141 

i  -  Sensor  Control  Unit 


1TBF 


Launches  Attempted 


Total  Aborts 


APPENDIX  D 


Computer  Proorams  for  Section  IV 


The  first  program,  ESTM5,  calculates  estimates  for  the  parameters 
mu,  sigma,  and  m  to  be  output  for  later  use  by  program  NEWTON.  This  pro¬ 
gram  uses  a  conjugate  gradient  search  algorithm  to  determine  the  zeros 
of  the  derivative  vectors  of  the  log  likelihood  equations  given  in  Sec¬ 
tion  IV.  The  actual  algorithm  is  the  IMSL  subroutine  2XCGR. 

The  input  should  be  on  TAPE  9  in  the  following  form:  dose,  number 
of  subjects,  and  number  of  positive  responses.  The  very  first  line  im¬ 
age  on  TAPE  9  should  only  contain  the  number  of  dose  levels.  The  user  is 
then  asked  to  input  initial  estimates  for  mu,  sigma,  and  m  by  the  inter¬ 
active  system.  The  output  is  written  on  TAPE  7  for  direct  use  by  the 
program  NEWTON,  or  prel iminary  inspection. 

The  second  code  listing  is  for  the  program  NEWTON.  This  program 
uses  a  Newton-Raphson  method  for  finding  the  zeros  of  the  derivative 
vectors  of  the  log  likelihood  functions.  The  program  expects  the  input 
data  to  be  on  TAPE  7.  The  input  data  should  be  in  the  following  order: 


1. )  First  card  -  N,  the  number  of  dose  levels 

2. )  Second  through  N+lst  card  -  dose,  number  of  subjects,  number  of 

responses 

3. )  N+2nd  card  -  initial  estimates  of  mu,  sigma,  and  m. 


The  final  estimates  of  mu,  sigma,  and  m  along  with  the  value  of  the 


log  likelihood  function  at  maximum  is  written  to  TAPE  8.  Also  on  TAPE  8 


is  the  information  matrix  evaluated  at  the  final  estimates. 


1*  PROGRAM  ESTM5  ***  THE  CONJUGATE  GRADIENT  SEARCH  METHOD*** 
2=  EXTERNAL  FUNLH 
3*  REAL  X(3),G<3),FLH,W<13) 

4=  COMMON/DA  TA/A(20,3MNUM,Z(20),DIFF<20) 

5*  READ(9,'(I3)')INUM 
6=  DO  222  IN*i,INUM 
7*  READ(9,*)A(IN,1)»A(IN,2)»A(IN,3) 

Z(IN)=ALOG(A(IN,l)) 

9*  DIFF(IN)*A(IN,2)-A(IN»3) 

10*222  CONTINUE 
11*  ACC*  .0000001 
12*  DO  111  1*1,3 
13*  PRINT*,'?' 

14*  READ*, X  (I) 

15=111  CONTINUE 

16*  CALL  ZXCGR(FUNLH»3»ACC,500»l,X,G,FLH,WtlER) 

17=  CALL  FUNLH(IDUM,X,FDUM,GDUM, -999 J 

18=  PRINT*,  FLH,X<1),X(2),X(3) 

19*  PRINT*, G(1),G(2),G<3> 

20*  END 

21*  SUBROUTINE  FUNLH(N»X,F,S,FLAG) 

22*  REAL  X(3),S(3),M2,F,P(20),DMU(20),DSIGMA  (20), DM1120) 

23*  COMMON/DA  TA/A(20,3)»INUM  ,2  (20),DIFF(20) 

24*  M2*  1.0 

25*  IF1FLAGJIE.-999JTHEN 
26*5  G1=GAMMA(X(3)) 


27*  G2*GAMMA(M2> 

28*  G3*GAMMA(X  (3)+M2) 

29*  BETA-G1*G2/G3 

^hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh. 

31»ftf**  THIS  DO  LOOP  CALCULATES  THE  CDF  (P1J))  FOR  EACH  DATA  POINT 
gg?«***  AND  THE  DERIVATIVE  HRT  EACH  OF  THE  PARAMETERS,  MU.SIGMA,  M  1 
gjffitt*  THE  CDF  FUNCTION 


35* 

36* 

37* 

38* 

39* 


41* 

42* 

43* 


45* 

46* 

47*10 

48* 

49* 

50*20 


HH 

DO  10  J*1,INUM 
Y*(Z(J)-X(i))/X(2> 

TEMP*EXP(Y)/(1+EXP(Y)) 

IF(X(3).LT.  0.0)  THEN 
P(  J)*  .9999999999 
ELSE 

P(J)»TEMP#*X<3) 

END  IF 

DM1  (J)«P(J)*ALOG(TEMP) 

PDF*  (EX  P(Y  #X  <3»*(  1 +EX  P<7  ))**(-!  (3)-M2)) /BETA 
DMU(J)*PDF/(-X(2)> 

DSIGMAl  J)*Y  *PDF/(-X  (2)) 

CONTINUE 
DO  20  1*1,3 
S(I)*0.0 
CONTINUE 


Mz*****«  ********** 4* 44 *4 *4 4**4 *4 *4 *4 44 *444444444444*  *4 *4*********** 
******* 

Elffi?*  rais  L00p  CALCULATES  the  derivative  VECTOR  S(I)  of  the  log- 

53- ******************** ************************ ********************* 
******* 

54*  DO  44  K-l.INUM 

55*  TEMP  1  *(DIFF  (K)-A  (K»2)*POD)  /  (P(D  *(i  -P(IC))) 

54*  S<l)«S<l)*TEMPl*DMUaO 

57*  S<2)*S(2)+TEMPl*DSIGMAOO 

58*  S(3)*S(3HTEMPi*DM10C) 

59*48  CONTINUE 
48*  FUNLL*8.8 

41*  DO  48  I-MNUM 

^kLL*FUNLL+DIFF(I)*ALOG(PU)/<l-P(I)))+A<I,2)*ALOG(l  -P(I» 

43*48  CONTINUE 
44*  F— FUNLL 

45*  S(l)*-S(l) 

44*  S(2)*-S(2) 

47*  S(3)=-S(3) 

48-  RETURN 
49*  ELSE 

78-  WRITEtt/dSniNUM 

71*  DO  188  H*i,INDM 

72-  HRITE(7>*)ZGCK)iA  GCZ«2)  »A  (IX  ,3) 

73*188  CONTINUE 

74-  WRITE(7,*)XU),X<2>,X<3> 

75*  END  IF 
74-  RETURN 
77*  END 


1*  PROGRAM  NEWTON 

2=  REAL  A(20,3>,P<20),DMU<20>,DS1GMA(20),DM1(20>,S(3),ZINFO<3,  3) 
3=  +X INV  (3 ,3) ,  WK(6 )  ,X  (3 )  >M2  >SS(3) 

4*  READ(7  »'(I3)  ONDATA 

5s  DO  100  IN=1,NDATA 
6=  READ(7»*)A(IN,l)fA(IN»2),A(IN,3) 

7*100  CONTINUE 
8*  READ(?,*)X(1),X(2),X(3> 

9*  EPSLON*  .000001 

10*  M2*  1.0 

11*5  G1*GAMMA(X(3)) 

12*  G2=GAMMA(M2) 

13*  G3*GAMMA  (X  (3)+M2) 

14*  BETA*G1«G2/G3 
15*  DO  10  J*i»NDATA 
16*  Y*(A(J»1)-X(1))/X(2) 

17=  TEMP*EXP<Y)/(1+EXP<7» 

18*  IPttO)  XT.  0.0)  THEN 

19*  P(J)=.9999999999 

20*  ELSE 

21*  P(J)=TEMP*#X  (3) 

22=  END  IF 

23*  DMl<J)=P(J)*ALOG(TEMP) 

24*  PDF*<EXP<Y#IQ))*(1+EXP(Y))**(-IQ)-M2))/BETA 

25*  DMU(J)*PDF/(-I(2)> 

26*  DSIGMA(J)*Y*PDF/(-X(2» 

27*10  CONTINUE 
28*  DO  20  1*1>3 
29*  S(I)*0.0 

30*  DO  30  J»  1,3 

31*  XINFO(I»J)*9.0 

32-30  CONTINUE 
33*20  CONTINUE 
34*  DO  40  X*1  iNDATA 

35*  TEMPl»((A(K,2)-A<r^))-A(K,2)*POD)/(PaD*(l-P(K))) 

36*  S(l)*S(l)+TEMPl*DMUOD 

37*  S<2)»S<2)+TEMP1*DSIGMA<K> 

38*  S<3)*S(3)+TEMP1  *DM1  (K) 

39*  TEMP2*AOC,2)/(PaD*(l-PaD» 

40*  XINPO(l,l>*XINFOa,l>+TEMP2*(DMUaC)**2) 

41*  X  INFO<  1  »2)*X  INFO  ( 1  *2)+TEMP2*«DMU  (K)*DSIGMA  GO)) 

42*  X  INFO(  1  »3)*X  INFO(l  »3)+TEMP2*((DMU  (K)*DM  1  GO)) 

43*  X  INFO  (2  »2)*X  INFO(2  >2)+TEMP2*(D  SIGMA  <D**2) 

44*  X  INFO(2»3)*X  INFO<2  »3)+TEMP2*«DSIGMA  G0*DM1  GD)) 

45*  XINFO(3^)»IINFO(3,3)+TEMP2*(DMl(r)**2) 

46*40  CONTINUE 

47*  XINFO(2,l)*XINFO(l,2) 

48*  XINFO(3»l)*XINFO(l»3) 

49*  X  INFO(3  »2)*I  INFO(2  >3) 

50*  TEST«SQRT(SU)*«2+S<2)**2+S(3)**2) 

51*  IF(TEST.GT.EPSLON)  THEN 

52*  CALL  LGINFCX INFO»3»3,3,0.0»X INV t3,SS,WlC,IER) 


DO  50  1*1,3 
Tl*Ti+XINV<l,I)*SU) 

T2*T2+X  INV(2,I)*S(I) 

T3*T3+XINV(3»I)*S<I) 

CONTINUE 
XCi)*X<i)+Ti 
X(2)*X(2)+T2 
X(3)*X<3)+T3 
GOTO  5 
ELSE 

PUNLL*0.0 
DO  60  I*1»NDATA 
ZZ*AU,2)-A(1,3) 

FUNLL*FUNLL+ZZ*ALOG(P(I)/(l-P(I)))+A(I,2)*ALQG(l-P(I) ) 
CONTINUE 

WRITE®, *)THE  LOG-LIKE  FUN  AT  MAX*  ',FUNLL 
WRITE®, *)'EST  OF  MU*  ',X<1) 

WRITE®, *)'EST  OF  SIGMA*  ',X<2) 

WRITE®»*)'EST  OF  Mi*  ',X<3) 

WRITE®, T  THE  FOLLOWING  IS  THE  INFORMATION  MATRIX")') 
DO  80  1*1,3 

WRITE®  ,*)X  IN  V  (1 , 1),X  INV(1 ,2)  ,X  INV  (1 ,3) 

CONTINUE 
END  IF 


appendix  e 


Computer  Programs  for  Section  V 

This  program  is  for  performing  the  calculations  of  the  Quantit  model 
of  Section  V.  This  code  is  almost  exactly  as  given  to  me  by  Mr.  Copen- 
haver  (Ref  11),  only  a  few  changes  have  been  made.  The  changes  made 
were  those  needed  in  order  to  implement  it  on  the  Cyber  759  computer  at 
WPAFB,  and  one  small  change  (lines  776-7 81)  to  handle  data  larger  than 
anticipated  by  Mr.  Copenhaver.  Copies  of  the  user's  manual,  also  ob¬ 
tained  from  Mr.  Copenhaver,  are  obtainable  from  Dr.  Barr  in  the  Depart¬ 
ment  of  Mathematics,  AFIT. 


74 


1=  PROGRAM  QUNTIT 
2=C**  QUANTIT  ANALYSIS 
3=C 

4=C***THIS  IS  THE  SECOND  FORTRAN  VERSION  OF  QUANTIT  ANALYSIS,  WRITTEN  FOR 
5=C  THE  IBM  370/159.  THIS  PROGRAM  WAS  CONVERTED  FROM  THE  PL/I  VERSION. 

6=C 

7=C***RESULTS  FOR  THE  NORMAL  MODEL  (PROBIT  ANALYSIS)  ARE  ALSO  PRODUCED  BY 
8=C  THIS  ROUTINE.  BOTH  THE  OMEGA  MODEL  (QUANTIT  ANALYSIS)  AND  THE  NORMAL 
9*C  MODEL  HAVE  BEEN  SCALED  SO  THAT  THE  PROBABILITY  DENSITY  FUNCTION  Fa) 
10=C  ATTAINS  A  MAXIMUM  VALUE  OF  1  AT  THE  ORIGIN  (F(0)  =  1). 
li=C  I.E., 

12=C  THE  OMEGA  PDF  IS:  Fa)  =  i  -  W«(V+1)  WHERE  W  IS  THE  ABSOLUTE  VALUE 
13=C  OF  (2*P-1). 

14=C  THE  NORMAL  PDF  IS:  F(X)  =  EXP(-Y*Y*PI)  WHIRE  PI=3.14159265 
15=C 

16=C*#THE  FOLLOWING  STATEMENT  DOUBLE  PRECISIONS  EVERY  VARIABLE  BEGIN  NING 
17=C  WITH  THE  LETTERS  A  THRU  H  AND  0  THRU  Z. 

18=  IMPLICIT  DOUBLE  FRECISION(A-H,0-Z) 

19=  DIMENSION  WHAT(7), XXX (50) 

20=  COMMON  /BLANX1/  XX(50),XN(50),XS(50),XP(50),XPHAT(50),XSAVE  P(50) 

21=  COMMON  /BLANK2/  AHAT,BHAT,AINIT, BINIT, V 

22=  COMMON  /BLANK3/  SAVEA(101),SAVEB(18i) 

23=  COMMON  /BLANK4/  ISETS,ITERA,ISWANA,IVOP,ISTEP,KD,MN,IXSW 

24=  COMMON  /BLANKS/  SAVEV(50),V123(3),SAVELN(50),XL123(3) 

25=  COMMON  /COM1/  XK,ID,XI 

26=  COMMON  /COM2/  PI,SQ2PI,A1,A2,A4,XLIKE 

27=  COMMON /ALFHAT/  VFIN,TITLE(9) 

28=  COMMON  /PANDV/  XPVAL(50),XVVAL(7) 

29=  CHARACTER  XCNTRL*5,XTITLE#5»XEDVAL*5,XDOSES*5,XVPARM*5,XFIN  IS#5, 

30=  1XBLANK#8,FINAL1*8,VFIN*5,TITLE*8,XLABEL*5 

31=  DATA  XCNTRL,XTITLE,XEDVAL,XD05ES,XVPARM,XFINIS  /'CNTRL',TI  TLE', 

32=  i'EDVAL', 'DOSES', 'VPARM'jTINlSV 

33=  DATA  FINAL1/  '  (FINAL)'/ 

34*  DATA  XBLANK  /'  7 

35=  PI  =  3.14159265359 
36=  SQ2PI  =  DSQRT(2.*PI) 

37=  ISETS  =  0 

38=C**  READ  CNTRL  CARD 

39=  1  READ(5,2,END*500)  ILABEL,KD,LOGT,ILOGA,MN,IVOP,INOV,IPRTV 

40=  2  FORMAT(A5,I2,I1,F3.0,I1,I1,I1,I1> 

41=  ISETS  =  ISETS  +  1 

42=  NODOS  =  0 

43=  IERR  =  0 

44=  WRITE(7»3)  ISETS 

45=  3  F0RMAT(1H1,'DATA  SET  M3) 

46=  IF  (XLABEL.EQJCCNTRL)  GO  TO  4 
47=  WRITE  (7^.) 

48=  5  FORMAT(1H0,'»  *  *  ERROR,  CNTRL  CARD  IS  NOT  PRESENT  OR  IS  NO  T  THE  F 

49=  1IRST  CARD  IN  THIS  DATA  SET.  PROGRAM  IS  TERMINATED') 

50=  GO  TO  500 

51*  4  CONTINUE 

52*  IFO.OGT.NE.3)  GO  TO  6 

53=C**CHECK  FOR  VALID  LOG  TRANSFORMATION 


54=  IF(XLOGA.GT.i  .0)  GO  TO  6 
55=  WRITE(7,7)  XLOGA 

56=  7  FORMATUH0,'*  *  *  ERROR  —  A=  ',F8.3,'  IS  INVALID  BASE  FOR  LOG  TRA 

57=  INFORMATION') 

58=  6  CONTINUE 

59=  IXSW  =  1 
60=  ISWANA  =  0 

61=  IEDCRD  =  0 

62=  IVCARD  *  0 

63=  VFIN  =  IBLANK 
64=C**READ  TITLE  CARD 

65=  READ(5»8 *END=500)  XLABEL,<TITLE(I>»I=i,9) 

66=  8  FORMA  T<A5,9A8) 

67=  IF(X LABEL.EQ JCTITLD  GO  TO  9 

68=  WRITE(7»10)  X  LABEL, (TITLE(I)»I= 1  ,9) 

69=  10  FORMAT(1H0,'*  *  *  ERROR:  TITLE  CARD  NOT  ENCOUNTERED.  THE  FO  LLOWING 
70=  1CARD  WAS  READ  IN:'»/iX»A5,9A8) 

71=  GO  TO  500 

72=  9  WRITE<7,1 1)  (TITLEU),1=1 ,9) 

73=  11  FORMAT(1H0,9A8) 

74=C«READ  NEXT  CARD 

75=  50  READ(5,12)  XLABEL»(WHAT(I)»I=1,7) 

76=  12  FORMAT(A5,5X,7Fi0.0) 

77=C**CHECK  IF  EDVAL  CARD 

78=  IF(XLABEL.NE.X  EDVAL)  GO  TO  13 

79=  IEDCRD  =  1 

80=  IF(  (MN.GE.1).AND.<MN.LE.7»  GO  TO  14 
81=  WR1TE17.15)  MH 

82=  15  FORMAT(1H0,'*  #  *  ERROR:  DATA  FROM  EDVAL  CARD  CANNOT  BE  RETR  IEVED. 
83=  1  CNTRL  CARD  INDICATES  THAT  THERE  ARE  ',14, 'VALUES') 

84=  IFUERR.EQ.0)  IERR  =  1 
85=  GO  TO  1000 
86=  14  DO  16  1=1, MN 
87=  16  XPVAL(I)  =  WHAT(I) 

88=  DO  17  1=  1»MN 

89=  IF((XPVAL(I).GT.0.)J1ND.(XFVAL(I).LT.1  J)  GO  TO  17 

90=  WRITE(7,18)  XPVAL(I) 

91=  18  FORMAT(1H0,'*  *  *ERROR:  THE  ED  VALUE  OF  P=  '.D12.5,'  IS  OUT  OF  RAN 
92=  1GE.  MUST  BE  BETWEEN  0  AND  1 ') 

93=  GO  TO  1000 

94*  17  CONTINUE 

95=C**CHECK  IF  VPARM  CARD 

96=  13  IFQC LABEL JIE.XVPARM)  GO  TO  19 

97=  IVCARD  *  1 

98=  IF<  (INOV.GE.l)^ND.(INOVlE.7) )  GO  TO  20 
99=  WRITE(7^1)  INOV 

100*  21  FQRMAT(1H0»'*  *  *ERROR:  DATA  FROM  VPARM  CARD  CANNOT  BE  RETR  IEVED 
101=  1COLUMN  13  OF  CNTRL  CARD  INDICATES  THERE  ARE  ',15,'  VALUES') 

102=  1FUERR.EQ.0)  IERR  =  1 
103=  GO  TO  1000 
104=  20  DO  22  1=1, INOV 
105=  XVVALU)  =  WHAT(I) 


106=  IF(  (XVVAL(I).GT.(-1.0)).AND.(XVVAL(I).LE.20.))GO  TO  23 
107=  WRITE<7,24)  XVVALU) 

108=  24  FORMA TUH0/*  *  *  ERROR:  INVALID  VALUE  OF  V=  SD12.5/  WAS  READ  FR 
109=  10M  VPARM  CARD.  MUST  BE  GREATER  THAN  -1  AND  LESS  OR  EQUAL  TO  20.') 

110=  GO  TO  1000 

lii=C**ONLY  FIRST  3  DECIMAL  PLACES  OF  V  ARE  USED. 

112=C***THE  FOLLOWING  3  STATEMENTS  ARE  IDENTICAL  TO  FLOORIX)  IN  PL/I. 

113*C  THAT  IS,  THE  LARGEST  INTEGER  .LE.TOX.  THE  FORTRAN  FUNCTIO  N 
114=C  IDINT(X)  IS  IDENTICAL  TO  FLOORtX)  FOR  X.GT.0  » BUT  NOT  FOR 
115=C  NEGATIVE  VALUES  OF  X. 

116=  23  CONTINUE 
117=  VDEL  =  IVVAL(I) 

1 18=  DELSGN  =  DSIGN(0.5D0,VDEL) 

119=  V  =  DBLE(  IDINT(VDEL*1000.  *  DELSGN)  )/l 000. 

120=  IF  (XWAL(I).LT.-.999)  XVVAUI)  =  -.999 
121=  22  CONTINUE 

122=C»CHECIC  IF  DATA  (l.E.  'DOSES')  CARD  HAS  BEEN  READ. 

123=  19  IF  (XLABEL.NE.X  DOSES)  GO  TO  25 
124=  NODOS  =  NODOS  +  1 
125=  IF  (NODOSXEJCD)  GO  TO  26 
126=  WR1TEC7.27)  XU 

127=  27  FORMAT(lH0/»  *  *ERROR:  THE  NO.  OF  DATA  CARDS  EXCEEDS  THE  V  ALUE  OF 
128=  1  ',13/  SPECIFIED  IN  COLUMNS  6-7  OF  CNTRL  CARD') 

129=  GO  TO  1000 

130=  26  XX (NODOS)  =  WHAT(l) 

131=  XS(NODOS)  *  WHAT12) 

132=  XN<NODOS)  =  WHAT13) 

133=  IF<  OCS(NODOS).LT.0.0).OR.GCS(NODOS).GTJ[N(NODOS)).OR.(XN(N  ODOS).L 
134=  1E.0.0))  GO  TO  29 

135=  GO  TO  30 

136=  29  WRITE!?, 31)  XX(NODOS),XS<NODOS),XN<NODOS) 

137=  31  FORMATUH0,'**  *ERROR:  ONE  OR  MORE  INVALID  DATA  ITEMS:  DOS  E=  ',D1 
138=  12.5/  S=  '.D12.5,'  N=  ',D12.5) 

139=  IF  (IERR.EQ.0)  IERR=1 

140=  GO  TO  1000 

141=  30  CONTINUE 

142=  25  1FQC LABEL .NE.XFINIS)  GO  TO  50 

143=  IF(  (MN.NE.0).AND.(IEDCRD.EQ.0))GO  TO  32 

144=  GO  TO  33 

145*  32  WRITE<7^4)  MN 

146=  34  FORMAT(1H0,'»  *  *ERROR:  MN=  ',15/  IN  COL.  12  OF  CNTRL  CARD  ,  BUT  N 
147=  10  EDVAL  CARD  IS  PRESENT.  DEFAULT  ED  VALUES  WILL  BE  USED') 

148=  GO  70  35 

149=  33  IF(MN.NE.0)  GO  TO  36 

150=C«*ASSIGN  DEFAULT  EDVALUES  IF  EDVAL  CARD  NOT  PRESENT 

151=  35  MN=7 

152=  XPVAL(l)  =.01 

153=  XPVAL(2)  =.05 

154=  XPVAL(3)  =.10 

155=  XPVAL(4)  =.50 

156=  IPVAL<5)  =.90 

157=  XPVAL(6)  =.95 


s 

I 


158*  XPVAL(7)  3 .99 

159*C*+CHECK  ON  OPTION  FOR  V  (COL.  13  ON  CNTRL  CARD) 

16«=C**CHECE  FIRST  FOR  SEARCH  PROCEDURE 
161s  36  IFUVOP.EQ.2)  GO  TO  37 
162=  IF(IVOP.EQ.e>  GO  TO  41 

163*C**CHECK  IF  VALUES  OF  V  ARE  READ  IN  (OPTION  2) 

164*  IF<  (IVOP.EQ.l).AND.(IVCARD.EQ.8)  )  GO  TO  38 

165*  GO  TO  39 

166*  38  WRITE(7»48)  IVOP 

167*  48  F0RMAT(1H8,'*  *  *ERROR:  IVOP  *  ',15,'  IN  COL.  13  OF  CNTRL  C  ARD,  BU 
168*  IT  NO  VPARM  CARD  IS  PRESENT.  DEFAULT  VALUES  ARE  ASSIGNED') 

169*  GO  TO  41 

178*  39  IF(IVOPJIE.e)  GO  TO  37 

1?1*C»ASSIGN  DEFAULT  VALUES  OF  V. 

172*  41  XVVAL(l)*  -.9 
173*  XWAL(2)*  -J5 

174*  XWALI3I*  8. 

175*  IWAL(4)*  1. 

176*  XWAL(5)*  5. 

177*  XWAL(6)*  18. 

178*  XWAL(7)*  28. 

179*  INOV  *  7 

188*  37  IF  (NODOS.Ea.ID)  GO  TO  1888 
181*  WRITE(7»43)  XD.NODOS 

182*  43  F0RMAT(1H8,'*  *  *ERROR:  ',15, 'DOSES  ARE  INDICATED  ON  CNTRL  CARD  (CO 
183=  1L.  6-7).  ONLY  ',15,'  DOSE  CARDS  WERE  PRESENT*) 

184*  IFUERR.EQ.8)  I  ERR  *  1 

185*  GO  TO  1888 

186*C«CHECK  IF  ERROR.  IF  SO, GO  TO  NEW  DATA  SET. 

187*  1888  IFUERR.EQ.1)  GO  TO  1 
188*  IF  (LOGT.EQ.4)  WRITE(7,44) 

189*  44  F0RMAT(1H8, DOSAGE  TRANSFORMATION:  NONE*) 

198*  1F(  (LOGTJLT.2).OR.(LOGT.GT.4) )  LOGT*  1 

191*  IF  (LQGT.EQ.1)  WRITE(7,45) 

192*  45  F0RMAT(1H8, DOSAGE  TRANSFORMATION:  LOG(BASE  18)') 

193*  IF  (LOGT.EQ.2)  WRITE(7,46) 

194*  46  F0RMAT(1H8, DOSAGE  TRANSFORMATION:  NATURAL  LN,BASE  EO 
195*  IFO.OGT.EQ.3)  WRITE(7,47)  XLOGA 

196*  47  F0RMAT(1H8, DOSAGE  TRANSFORMATION:  LOG(BASE  ',F5.8,' )') 

197*  WRITE(7,48) 

198*  48  F0RMAT(1H8,17X, 'TRANSFORMED  NO.  OF  NUMBER*) 

199*  WRITE(7,49) 

288*  49  FORMATdH  ,7X , DOSAGE', 7X, 'DOSE', 8X, 'SUBJECTS' ,4X,'RESPONDI  NG',4X, 
281*  1 'PROPORTION*) 

282*  1*8 
283*  162  1*1+1 

284*  XSAVEP(I)  *  XS(I)/XN(I) 

285*  XXX(l)  *  XX(I) 

286*  XP(I)  *  XSAVEP(I) 

287*  IF(  (XS(I).EQ.8).OR.(XS(I).EQ.XN(I)) )  XP(I)  *  (XS(I)+  8^)  / 

288*  lGCN(I)  +  1.8) 

289*  IFO.OGT.EQ.4)  GO  TO  8881 

218*  8888  IF  (  LOGT.EQ.1)  XKI)  *  DLOG18(XX(I)) 
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f  f 


211*  1F(LQGT.EQ.2)  XX(I)  *  DLOGUCXd)) 

212*  IF(LQGT.NE.3)  XLQGA  *  13. 

213*  IF(LOGT.NE.3)  GO  TO  3001 

214*  XTRANS  *  DLOG10(XLOGA) 

215*  XX(1)  *  DLOG10GCX  (I)) /XTRANS 

214*  8001  CONTINUE 

217*  WRITEC7»43)  XXX(I),XX(I),XN(I),XS(I>»XSAVEP(I) 

218*  63  FORMAT(1H0,3X,2(D12.5,2X)»F6.0»7X,F6.0,9X»F9.6) 

219*  IFd.LT.XD)  GO  TO  142 

220*  WRITE(7,51>  XD 

221*  51  FORMA  T(1H0//,'  NO.  OF  DOSAGE  LEVELS  *  ',I4> 

222*  WRITE(7,52) 

223*  52  FORMAT(1H0,'THE  OPTION  FOR  V  *  ') 

224*  IFdVOP.NE.2)  GO  TO  53 
225*  WRITE(7,54> 

224*  54  FORMAT!  1H+»20X /SEARCH  PROCEDURES 
227*  IFdPRTV.EQ.1)  GO  TO  55 
228*  WRITE(7,54) 

229*  54  F0RMAT(1H+»37X, '(PRINT  COMPLETE  RESULTS  FOR  V=0,1,  AND  FINA  L  V)') 
230*  GO  TO  60 

231*  55  WRITE(7,57) 

232*  57  FORMA  T(1H+»37X /(PRINT  COMPLETE  RESULTS  FOR  ALL  VALUES)' ) 

233*  GO  TO  60 

234*  53  CONTINUE 

235*  IFGVQP.EQ.1)  GO  TO  58 

236*  WRITE(7»59)  (XWAL(I),I*i,INOV) 

237*  59  FORMAT(1H+*20X »'  DEFAULT  VALUES:  V*  ',  7(F10.3) ) 

238=  GO  TO  60 

239=  58  HRITE(7f41)  (IVVAL<I>,I*l,INOV) 

240*  61  FORMAT(1H+»20X ,'  INPUTTED  VALUES:  V*  ',  7(F10.3) ) 

241*  60  CALL  MLEAB 

242*  CALL  PRINT 

243*  ISWANA  *  1 

244*  IFUVOP.EQ.2)  GO  TO  65 

245*C»PRODUCE  RESULTS  FOR  FIXED  VALUES  OF  V  (IVOP  *  1  OR  2) 

246*  DO  66  I*l,INOV 
247*  V*  XVVAL(I) 

248*  CALL  VRAT 

249*  CALL  MLEAB 

250*  CALL  PRINT 

251*  66  CONTINUE 
252*  WRITE(7,3000)  ISETS 

253*  3000  FQRMAT(1H0»'*  *  *  E  N  D  OF  DATA  SE  T',I3,'  *  *  *0 
254«C**GO  TO  NEW  SET  OF  DATA 
255*  GO  TO  1 

256eC** SEARCH  PROCEDURE:  FIND  V  IN  (-1  <  V  <*  20  )  THAT  MAXIMIZES  THE  LIKELI 

257*  65  ISTEP  *  0 

258*  IVDONE  *  3 

259*  DO  67  1*1,3 

260*  ISTEP  *  ISTEP  +  1 

261*  TEMPV  *  DBLE(I)  -  1. 

262*  V  *  TEMPV 
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263*  V123(I)  *  TEMPV 

264=  SAVEV(I)  =  TEMPV 

265=  CALL  VRAT 

266=  CALL  MLEAB 

267=  XL123G)  =  XLIKE 

268=  SAVELNd)*  XLIKE 

269=  IP  ( (I.LT.3).0R.(IPRTV.EQ.i) )  CALL  PRINT 

270=  IF  (  d.EQ.3).AND.dFRTV.EQ.0))GO  TO  68  271=  GO  TO  67 

272=  68  WRITE(7, 70) 

273=  70  FORMATdH  1  »58X /SEARCH  PROCEDURE', //IX, 12X, 'STEP', 12X,,V,,1  2X,TN 
274=  1L' ) 

275=  WRITE<7,7?)  ISTEP,V, XLIKE 

276=  67  CONTINUE 

277=  IF<  (XL123(2).GT.XL123(1)).AND.(XL123(2).GT.XL123(3))  )  GO  TO  72 

278=  GO  TO  73 

279=  72  DEL  =  -0.4 

280=  DO  74  I  =  1,2 

281=  DEL  =DEL/2. 

282=  DO  75  J=l,2 
283=  ISTEP  =  ISTEP  ♦  1 
284=  IVDONE  =  IVDONE  +  1 
285=  VDEL  =  V123(2)  ♦  DEL 

286=  DELSGN  =  DSIGN<0.5D0,VDEL) 

287=  V  =  DBLEl  IDINT<VDEL*1000.  +  DELSGN) )/ 1000. 

288=  V123(2*J-1)  =  V 

289=  SAVEVdVDOND  =  V 

290=  CALL  VRAT 

291=  CALL  MLEAB 

292=  XL123(2*J-1)  =  XLIKE 

293=  SAVELNdVDONE) »  XLIKE 

294=  IF(IPRTV.EQ.l)  GO  TO  175 

295=  WRITE(7,77)  ISTEP ,V,X  LIKE 

296*  77  FORMAT(1H0,11X»I4,10X,F7 .3,4X4)15.7) 

297=C**FORMAT  77  IS  FI  IN  PL/I 

298=  GO  TO  176 

299=  175  CALL  PRINT 

300*  176  CONTINUE 

301=  DEL  «  -DEL 

302=  75  CONTINUE 

303=C*GO  TO  HILQW 

304=  IF  ( (XL123(3).GT.XL123<2)).OR.(XL123(l).GT.XL123<2»  )  GO  TO  78 
305=  74  CONTINUE 
306=  V*l. 

307=C**GQ  TO  VI 

308=  GO  TO  80 

309=  73  CONTINUE 

310=C**THIS  IS  HILOW 

311=  78  IF  (XL123(3).LE.XL123<2))  GO  TO  81 

312=  lFdVDONEJlE.3)  GO  TO  82 

313=  DEL  =  2. 

314*  GO  TO  83 

315*  82  DEL* -DEL 

316=  83  XL123<1)  =  XL123<3) 


317=  V123(l)  =  V123(3) 

318=  GO  TO  84 

319=  81  IF  (IVDONE.EQ.3)  DEL  =  -.2 

320=C**DELV 

321=  84  CONTINUE 

322=  VDEL  =  V123(l)  ♦  DEL 

323=  DELSGN  >  DSIGN(0.5D0,VDEL) 

324=  V  *  DBLE(  IDINT(VDEL*1000.  +  DELSGN)  )/1000. 

325=  7123(2)*  V 

326=  ISTEP  =  ISTEP  +  1 

327=  DO  85  I  =  1, IVDONE 

328*  IF(V.NE.SAVEV(I))  GO  TO  85 

329*  IFUPRTV.EQ.0)  GO  TO  86 

330=  WRITE(7,87)  ISTEP^AVEV(I) 

331=  87  FORMAT(1H0,//,1X,'STEPM3,':  V=  ',F5.2,'  (PREVIOUSLY  CALCU  LATED)' 
332*  1) 

333=C**FORMAT  87  IS  F2  IN  PL/I 
334=  GO  TO  88 

335=  86  WRITE(7,77)  ISTEP,SAVEV(I),SAVELNG) 

336=  88  XL123(2)  =  SAVELN(I) 

337=C**GO  TO  CHKL 

338=  GO  TO  89 

339*  85  CONTINUE 

340=  IVDONE  =  IVDONE  +  1 

341*  SAVEVdVDONE)  *  V 

342=  CALL  VRAT 

343=  CALL  MLEAB 

344=  IFdPRTV.EQ.1)  GO  TO  91 

345=  WRITE(7,77)  ISTEP, V,XLIKE 

346=  GO  TO  92 

347*  91  CALL  PRINT 

348=  92  XL123(2)  =  XLIKE 

349=  SAVELNdVDOND  =  XLIKE 

350»C**CHKL 

351*  89  IF(XL123(2).GT.XL123(1))  GO  TO  93 
352*  DEL  »  -DEL/2. 

353*  XABS  *  DABSCDEL) 

354*  IF(  (XABS.GE.0.5)  .AND.(V.GT.5.0) )  GO  TO  94 

355*  IF«XABS.GE.0.25).AND.(V.GT.2.0).AND.<V.LE.5.0) )  GO  TO  94 

356*  IF(  <XABS.GE.0.1).AND.<V.LE.2.0) )  GO  TO  94 

357*  GO  TO  95 

358*  94  XH23(1)  »  XL123(2) 

359*  V123(1)  *  VI  23(2) 

360*C**GO  TO  DELV 
361*  GO  TO  84 

362*  95  CONTINUE 
363*  GO  TO  96 

364*  93  CONTINUE 

365*C«*IF  V*-,8  AND  DEL*-.2,  SET  DEL  *  -.1 
366*  IF<  DABS(V+.8D0).GT.i.0D-04)  GO  TO  8002 
367*  IF(  DABS(DEL+0.2D0)lE.1.0D-04  )  DEL  =  -0.1D0 
368*  8002  CONTINUE 


369*C***CHECK  IF  V*-.?  OR  V*20. 

370*  IF(  DABS(V+0.9D0).LE.1.8D-08)  GO  TO  100 
371*  IFCV.GE.20.D0)  GO  TO  100 
372=  XL123<1)=XL123<2) 

373*  V123<1)  *  V 

374»C**GO  TO  DELV 
375*  GO  TO  84 
376*  94  V*  V123(l) 

377»C**V1 

378*  80  CALL  VRAT 
379*  CALL  MLEAB 

380*  100  VFIN  *  FINAL1 
381»C**FINALV 
382*  CALL  PRINT 

383*  HRITE(7.3000)  ISETS 

384*  WRITEC7,101> 

385*  101  FORMAT(iHl) 

386=C»*GO  TO  NEW  DATA  SET 
387=  GO  TO  1 

388*  500  HR1TE(7|501) 

389*  501  FORMAT(1H0»//»20X»'*  #  *  E  N  D  OF  PROGRAM***') 

390*  STOP 

391*  END 

392*  SUBROUTINE  PRINT 

393*  IMPLICIT  DOUBLE  PRECIS10N(A-H,0-Z) 

394*  DIMENSION  PDIF(50),EDP(50)»VAREDP(50)»SEEDP(50) 

395*  COMMON  /PANDV/  XPVAL(50),XWAL(7> 

396*  COMMON  /BLANK1/  XX(50)  JN(50),XS<50),XP(50),XPHAT(50),XSAVE  P(50) 

397*  COMMON  /BLANX2/  AHAT»BHATfAINIT,BINIT»V 

398*  COMMON  /BLANK3/  SAVEA(101),SAVEB<i01) 

399*  COMMON  /BLANK4/  ISETS, ITERA,ISWANA,IVDP,ISTEP,XD,MN,IXSW 

400*  COMMON  /BLANK5/  SAVEV(50),V123(3)^AVELN(50),XL123(3) 

401*  COMMON  /ALPHAT/  VFIN,TITLE(9) 

402*  COMMON  /COM2/  PI,SQ2PI*AlfA2»A4 ,XLIKE 

403*  CHARACTER  XBLANX*8,VFIN*5,TITLE*8 

404*  DATA  XBLANK  /'  7 

405*  I  TENS*  -8 

406-  CHISQ  *  0. 

407*  IFUSHANA.EQ.9)  GO  TO  1 

408*  WRITE(7,2> 

409*  2  FORMATOHl »50X ,'*  *  *  QUANTIT  ANALYSIS  *  *  *') 

410*  IF(  (IVOPJEQ.2).AND.<VFIN.Ea.XBLANXD)GO  TO  3 
411*  WRITE(7»4> 

412*  4  FOHMAT(1H0) 

413*  GO  TO  5 

414*  3  WRITE  (7 ,6)  ISTEP 

415*  6  FORMAT(1H0, 'SEARCH  PROCEDURE:  STEFM3) 

416*  5  WRITE(7»7)  V,VFIN 

417*  7  FORMATdH  ,61X  ,*V  *  ',F7.3,A8) 

418*  GO  TO  8 

419*  1  WRITE(7»9) 

420*  9  F0RMAT(1H1,/51X,'*  »  *  PROBIT  ANALYSIS  *  *  *') 

421*  8  CONTINUE 


422=  WRITECM0)  1SETS,(TITLEU),I=1,9) 

423=  10  FORMATdH0,'DATA  SET  ',13,':  ',9A8  > 

424=  WR1TE(7,1 1)  A1NIT»B1NIT 

425=  11  FORMATdH0,/, IX, 'INITIAL  ALPHA  =  ',D15.7,7X, 'INITIAL  BETA  =  ',D15. 

42 6=  17) 

427=  WRITE<7,12> 

428=  12  PORMATdHO, 'ITERATION', 12X,'ALPHA',14X, 'BETA 0 

429=  DO  13  1=1,ITERA 

430=  WRITEC7.14)  I,SAVEA(I),SAVEB(I) 

431=  14  FORMATdH  ,2X,I3,UX,D15.7,3X,D15.7) 

432*  13  CONTINUE 

433=  AA2  >  -A  2 

434=  HR1TE(7,15)A4,AA2,AA2»A1 

435=  15  FORMATdH0,/iX,'THE  VARIANCE-COVARIANCE  MATRIX  FOR  ALPHA  AN  D  BETA 
436=  l',/,2(/2(6X,D15.7)) ) 

437=C**N_LINE 

438=  16  ITENS  =  ITENS  +9 

439=  NUM  =  MlN0OCD,ITENS+8) 

440=  WRITE(7,17)  <XPHAT(I)»I=ITENS,NUM) 

441=  17  FORMATdH0,//lX,'MLE  FOR  P:  ',9(3X,F10.7> ) 

442=  WRITE(7,18)aSAVEP(I),I=ITENS,NUM) 

443=  18  FORMATdH  , 'OBSERVED  P: ',9(3X,F10.7)  ) 

444=  DO  19  1=1  JO) 

445=  19  PDIF(I)  =  XPHAT(I)  -  XSAVEP(I) 

446=  WRITE(7,20>  ( PDIF(I),I=ITENS,NUM) 

447*  20  FORMATdH  .DIFFERENCE  s  ',9(3X,F10.7)  ) 

448=  IF(NUMJ.TJCD)  GO  TO  16 

449*  BHATSQ  =  BHAT*BHAT 

450=  IFIRST  =  A4  /  BHATSQ 

451=  ILAST  =  -A2*2.0  /  BHATSQ 

452=  XSEC  =  A1  /  BHATSQ 

453=  IF  (ISHANA.EQ.0)  GO  TO  21 

454=C«QUANTIT  EDVALUES 

455=  DO  22  I=1,MN 

456=  CALL  QUANT1  (IPVAL(I),EDP<I) ) 

457=  22  EDP(I)  =  (EDP(I)  -  AHAT)  /  BHAT 

458=  GO  TO  23 

459*C*»PROBIT  RESUTLS 

460=  21  CALL  INVNOR(X  PVAL,EDP,MN) 

461*  DO  24  I=1^1N 

462*  24  EDP(I)  *  <EDP<I)/  SQ2PI  -  AHAT)  /  BHAT 
463=  23  DO  25  1=1  tMN 

464=  VAREDP(I)  =  XFIRST  +  XSEC#EDP(I)*EDP(I)  +  EDP(I)*ILAST 
465=  25  SEEDP(I)  *  DSQRT  (  VAREDP(I) ) 

466=  WRITE<7,26) 

467*  26  FORMAT(1H0,/6X  ,'P',9X  ,'ED  ESTIMATE', 10X , 'VARIANCE', 10X  ,'STD  .  ERROR 
468=  19 

469*  DO  27  1=1, MN 

470=  27  HRITE(7,28)  XPVAL(I),EDP«),VAREDPa),SEEDP(I) 

471=  28  FORMA TdH0,F10.7^X,D15.7,2(4X,D15.7) ) 

472=  DO  29  1=1,0) 

473=  PDIF(I)  =  PDIF(I)*PD1F(I) 


474=  XYZ  =  XPHAT(I)*  <1.  -  XPHAT(I) ) 

475=  2?  CHISQ  =  CHISQ  +  XN(I)  *  PDIF(I)  /  XYZ 
476=  WRITE(7,30)  XLIKE, CHISQ 

477=  30  FQRMAT(1H0»/»5X,'LN  L  =  ',D15.7,  5X/CHI  SQUARE  =  ',D15.7  ) 

478=  RETURN 

479=  END 

480=  SUBROUTINE  QUANCD(XPHAT»XX  ,F> 

481=C**THIS  ROUTINE  CALCULATES  THE  CDF  FOR  THE  OMEGA  DISTRIBUTION.  I.  E., 
482=C  GIVEN  X,  FIND  P. 

483*  IMPLICIT  DOUBLE  PRECISION  (A-H,D-Z) 

484=  COMMON  /BLANK2/  AHATfBHAT*AINIT»BlNIT»V 

485=  COMMON  /BLANK4/  ISETS,ITERA,ISWANA,IVOP,ISTEP,ED,MN,IXSW 

486*  DIMENSION  IFHAT(50),XX(50>,F<50) 

487*  DO  1  I*  1»KD 

488*  2  *  AHAT  +  BHAT  *  XX(I) 

489=  K*0 

490*  Q2  =  1. 

491=  DO  2  J  =  1,101 

492=  PC  *XPHAT(I) 

493=C**ANOT 

494*  3  IF  (PCiT.0.5)  PC  =  1.0  -  PC 

495=  CALL  QUANT UPCfHP) 

496*  G  =  HP  -  DABS(Z) 

497=  Q1  *  1.  -  <DABS(2.*PC-1.))#*<V+1.) 

498*  XPHAT(I)  =  PC  -  G*Q1 
499=  IF  <XFHAT«).LT.1.0>  GO  TO  4 

500*C«FIND  NEW  INITIAL  ESTIMATE  OF  P.  IF  INITIAL  IS  GREATER  THE  FINA  L 

50i=C  ESTIMATE  ( FOR  P  GREATER  THAN  JS)  ,  THEN  CONVERGENCE  IS  GUARA  NTEED. 

502*  +  l 

503*  IF  OCEX/T.l)  GO  TO  5 

504*  PC  *  .9999D0 

505*C«K}O  TO  ANOT  AND  TRY  AGAIN 

506*  GO  TO  3 

507*  5  IF  (KK.GT.2)  GO  TO  6 

508*  PC  *  .99999999D0 

509*  GO  TO  3 

510=C»NOTE  THAT  CONVERGENCE  CRITERIA  FOR  P  IS  0.00001  .  HENCE  IF  T  HE 

511=C  PROGRAM  REACHES  THIS  POINT,  SET  P  =  .99999999 

512=  6  XPHAT(I)  =  .99999999D0 

513*C»G0  TO  FIN 

514*  GO  TO  10 

515*  4  Q2  *  DABSt  XPHAT(I)  -  PC  ) 

516*  IF(Q2XT.0.00001)  GO  TO  10 

517=  IF  (JJ.E.100)  GO  TO  2 
518*  WRITE(7,11)I 

519*  11  FORMAT(1H0,/NOTE:  MORE  THAN  :100  ITERATIONS  ARE  REQUIRED  FOR  P(', 
520=  113,0  IN  SUBROUTINE  QUANCD' ) 

521*  GO  TO  10 
522=  2  CONTINUE 

523=C«FIN 

524*  10  IF(Z.LT.0.0)  XPHAT(I)  *  1.  -  XPHAT(I) 

525=  1  F(I)  «  1.  -  (DABS(2.*XPHATU)  -  1.)  )**(V+1J 

526=  RETURN 
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527*  END 

523*  SUBROUTINE  NQRMCD  <XPHAT,XX,KD,F) 

529=C**TH1S  ROUTINE  CALCULATES  THE  CDF  FOR  THE  NORMAL  DISTRIBUTION.  I  .E., 
530*C  GIVEN  X,  FIND  P. 

531*  IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

532*  COMMON  /BLAN1C2/  AHA  T,BHAT,AINIT, BINIT,  V 

533*  COMMON  /COM2/  PI,SQ2PI,Al,A2,A4,XLirE 
534*  DIMENSION  XPHAT(50), XX  (50), F(50) 

535*  DIMENSION  B(5) 

534*  B0  *  0-2316419 

537*  B(l)  *  0.31938153 

538*  B(2)  *-0.356563782 

539*  B(3)  *  1.781477937 

540*  B(4)  *-1.821255978 

541*  B(5)  *  1.330274429 

542*  DO  1  1*1,  KD 

543*  CDF  =  0. 

544*  2  *  (AHAT  +  BHAT  *XX(I) )  *  SQ2P1 

545*  T  *  1.0/  (1.0  +  B0*DABS(D  ) 

546*  DO  2  J*  1,5 

547*  2  CDF  *  CDF  +  B(J)*T**DBLE(J) 

548*  CDF  *<CDF  /  SQ2P1)  *  DEXP(-Z*2/2J 
549*  XFHAT(I)  *  1.  -  CDF 

550*  IFG.LT.0.0)  XPHAT(I)  *  CDF 

551*  Z  *  AHAT  +  BHAT*XX(I) 

552*  1  F(I)  *  DEXP(-Z*Z*PI) 

553*  RETURN 

554*  END 

555*  SUBROUTINE  INVNOR  (P,Y,N) 

556*C**THIS  SUBROUTINE  CALCULATES  THE  NORMAL  DEVIATE  (I .E., MODIFIED  P  RDBIT) 
557*C  OF  P.  GIVEN  P,  FIND  T. 

558*  IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z> 

559*  DIMENSION  P(50),Y(50),CC(2),DD(3) 

560*  CC(1)  *  0.802853 

561*  CC(2)  *  0.010328 

562*  DD(1>  *  1.432788 

563*  DD<2)  *  0.189269 

564*  DD(3)  *  0.001308 

565*  DO  1  I  *  i,N 

566*  XNUM  *  2.515517 

567*  XDEN  *  1.0 

568*  PP  *  PCI) 

569*  IF(P(I).GT.0.5>  PP  *  1.  -  P(I) 

570*  T*  DLOG(l  .0  /  (PP*PP)  > 

571*  T*  DSQRT(T) 

572*  DO  2  J  *  1,2 

573*  2  XNUM*  XNUM  +  CC(J)*T**DBLE(J) 

574*  DO  3  J*  1,3 

575*  3  XDEN  *  XDEN  ♦  DD(J)  *  T»DBLE(J) 

576*  ZP  -  T  -  XNUM/IDEN 

577*  YU)  *  -ZP 

578*  IF  <PH).GT.0.5)  Y«)  *  ZP 
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579* 

1  CONTINUE 

l: 

■;! 

580* 

RETURN 

k- 

S  >v- 

581* 

END 

■r 

i  : 

582= 

SUBROUTINE  VRAT 

1 

583=C**THI5  ROUTINE  EXPRESSES  V+i  AS  A  RATIONAL  NUMBER.  THIS  IS  NECES  SARY 

584=C 

FOR  CALCULATION  OF  THE  QUANTIT  <IN  SUBROUTINE  QUANT1) 

s 

V- 

585* 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

r' 

V. 

586* 

COMMON  /COM1/  XK,XD,XI 

£ 

587= 

COMMON  /BLANK2/  AHAT,BHAT,AINIT, BINIT, V 

i 

588* 

DIMENSION  A<2,2) 

1 

1 

589»C**FIND  RATIONAL  NUMBER  F*V+i  =  XK/XD 

*,' 

590* 

F*  V+l 

m' 

591* 

W1  *  DBLEHDINTtF) ) 

V 

592= 

X*  F  -  W1 

L* 

\ 

593* 

X*  DBLE(  IDINT(X*1000.  +  .5) ) 

i 

% 

i 

594* 

Y=  1000. 

1 

I 

595* 

A(l,l)  =  Wi 

;< 

»,■ 

596= 

A  (1,2)=  1. 

r 
<  * 

597* 

A  (2,1)*  1. 

V 

598* 

A  (2,2)=  0. 

* 

V 

599=C**IF  X=0  ,  GO  TO  CALC 

he 

a 

600* 

IF(X.EQ.0.)  GO  TO  2 

fl 

p 

601* 

B  =  DBLE(  IDINT  (Y/X) ) 

- 

602* 

CALL  MULT(A,B) 

> 

603* 

INUM*  Y 

V 

604* 

DENOM  *  X 

.* 

605- 

DO  1  1  =  1,200 

a 

i  a 

•* 

606* 

WORK1  * INUM 

7 

607* 

XNUM  -  DENOM 

kT 

•«. 

608* 

DENOM  *  DMOD(WORKi, DENOM) 

v‘ 

IS 

609* 

IF  (DENOM.EQ.0.)  GO  TO  2 

V 

%' 

610* 

B*  DBLE(  IDINT(XNUM/DENOM) ) 

•r. 

611* 

1  CALL  MULT(A,B) 

1 

■ 

B 

612=C#*CALC 

613* 

2  XK  *  A(l,l) 

1  M 

* 

614* 

XD  *  A (2,1) 

5; 

-- 

615* 

XI  *  -1 

V 

616* 

TWO  *2. 

S 

*  B 

•* 

617* 

IF(DMOD<XE,TWO).EQ.0.)  XI  *  1 

i 

■ 

618* 

RETURN 

< 

\ 

619* 

END 

i 

620* 

SUBROUTINE  MULT(A,B) 

,T 

'! 

621* 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

■-* 

622* 

DIMENSION  A(2,2),C(2,2) 

t'- 

623* 

C(l,l)  *  A(1,1)*B  ♦  A(l, 2) 

fl 

i 

624* 

C(1 ,2)  *  A(l,l) 

-* 

625* 

C(2,l)  *  A(2,1)*B  +  A (2, 2) 

** 

> 

626* 

C(2,2)  *  A (2,1) 

>* 

627* 

DO  1  1*1,2 

•* 

628* 

DO  1  J*l,2 

A 

•* 

629* 

1  A(IJ)  *  C(I,J) 

n 

i 

630* 

RETURN 

i** 

r 

r 

■  ■ » 

*  . 

v> 

* 

631* 

END 

r 

* 

i* 

h 

</ 

■ 
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632=  SUBROUTINE  QUANT  1  (P.HP) 

633=C»*THIS  ROUTINE  CALCULATES  THE  QUANTIT  OF  P.  I  .E., GIVEN  P,FIND  H(  P). 
634=  IMPLICIT  DOUBLE  PRECISION  <A-H,D-Z> 

635=  COMMON  /COMi/  XK,XD,XI 

636=  COMMON  /COM2/  PI,SQ2PI,A1,A2,A4,XHKE 

637=  COMMON  /BLANK2/  AHAT,BHAT,AINIT,BINIT>V 
638=  P*V+1. 

639=  IF  (P.NE.8.5)  GO  TO  1 

648=  HP  =  8. 

641=  RETURN 

642=  1  PP  =  DABS(  2.*P  -1 J 

643=  PP1=  DEXP(  DLQG(PF)/XD) 

644*  PP2=  DEXP(  2.*DLOG<PP)/XD) 

645=  SUM1  =  8. 

646=  SUM2  =  8. 

647*  SUM3  =  8. 

648*  IMAX  =  IDINT<  OCK-l.)/2.) 

649=  IF  (IMAXJLT.l)  GO  TO  288 

658=  DO  28  1=1, IMAX 

651=  TEMPI  =  DBLE(I) 

652=  WORXi  =  DCOS(2.*PI*TEMPI*XD/XK) 

653=  WORK2  =  DCOS(2.*PI*TEMPI/XK) 

654=  WORK3  =  DSIN(2.*PI*TEMPI*XD/XK) 

655=  WORK4  =  DSIN(2.*PI*TEMPI/XK> 

656*  SUM1  =  SUM1  -  CID/XK)  *  WORK!  *  DLOG(l.-  2.*PPl*WORX2  +  PP2  ) 

657=  WORKS  =  PPUWORK4/C1.  -  PPl*WORK2) 

658*  28  SUM2  *  SUM2  +  2.*<XD/XK)#WORK3  *  DATAN(WORKS) 

659=  288  CONTINUE 

668=  IMAX  =  IDINT(  <XD-i.)/XK) 

661=  IF  (IMAX.LT.1)  GO  TO  281 

662'  DO  38  I  *  1,IMAX 

663=  TEMPI  =  -DBLE(I) 

664=  WORK1  =  TEMPI  *  F  +  1. 

665=  38  SUM3  =  SUM3  -  DEXPI  WORK1  *  DLOG(PP)  -  DLOG(WORKl> ) 

666=  281  HP  =  SUM1  +  SUM2  ♦  SUM3  -  (XD/XK)*DLOG(i.-PPl) 

667=  1  +  XD*(l.+XI)*DLOG(l.  +  PF1)/(2.*XD 

668=  HP  =  HP/2JD8 

669=  IF  <P.LT.8.5>  HP  =  -HP 

678*  RETURN 

671=  END 

672=  SUBROUTINE  MLEAB 

673=C**THIS  ROUTINE  CALCULATES  THE  MLE'S  FOR  ALPHA, BETA 
674*  IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

675=  COMMON  /BLANK1/  XXI58)  JN(58>,XS(5W,XP(58),XPHAT(58),XSAVE  F(58> 
676=  COMMON  /33LANK2/  AHAT3HAT^INIT,BINIT,V 

677=  COMMON  /BLANK3/  SAVEA(181),SAVEB(18i) 

678*  COMMON  /BLANK4/  ISETS,ITERA,ISWANA,IVOP,ISTEP,KD,MN,IXSW 

679=  COMMON  /BLANK5/  SAVEV(58),Vi23(3),SAVELN<58).XLi23C3> 

688*  COMMON  /COM2/  PI,SQ2PI,A1,A2,A4,XLI1CE 
681=  DIMENSION  7(58), F(58) 

682=  EPSI  =  8.881 

683=  MAX  IT  =  188 
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484=  HIKE  =  9. 

<585=  ITERA  =  0 

484=  IFUSWANA.EQ.0)  GO  TO  1 

487=  DO  2  1=1, KD 

488=  CALL  QUANT1(XP(I),Y(I» 

489=  2  XPHAT(I)  =  XP(I) 

490=  GO  TO  3 

491=  1  CALL  INVNOR(XF,Y  »KD) 

492=C« 

493=  DO  20  I  =  i,KD 
494=  20  Y  (I)  =  Y  (D/SQ2PI 

495=C**INITIAL  ALPHA  ,BETA  BY  LEAST  SQUARES 
494=  3  IF(IXSW.EQ.0)  GO  TO  4 

497=  SUMX  =  0. 

498=  XTX  =0. 

499=  DO  5  1  =  1  ,KD 

700=  SUMX  =  SUMX  +  IX  (I) 

701=  5  XTX  =  XTX  +  XX(I)*XXU) 

702=  SSX  =  XTX  -  (SUMX*SUMX)/DBLE£KD) 

703=  XBAR  =  SUMX/DBLEUCD) 

704=  IXSW  =  0 

705=  4  CONTINUE 

704=  XTY  =  0. 

707=  SUMY  =  0. 

708=  DO  4  1  =  1,KD 

709=  XTY  =  XTY  +  XX(I)*Y(I) 

710=  4  SUMY  =  SUMY  +  Y(I) 

711=  YBAR  =  SUMY/DBLEUCD) 

712=  SSP  =  XTY  -  (SUMX  *  SUMY)/DBLE(KD) 
713=  BHAT  =  SSP/SSX 

714=  AHAT  =  YBAR  -  BHAT*XBAR 

715=C**NEWTON-RAPHSON  PROCEDURE 
714*C 

717*C»ITER 
718*  W?ITE(7,51) 

719=  51  FORMAT  ( 1H0  ,'**  1  «*0 

720*  7  IFUTERA.EQ.0)  GO  TO  8 

721=  SAVEAUTERA)  =  AHAT 

722=  SAVEBUTERA)  =  BHAT 

723=  GO  TO  9 

724=  8  AINIT  =  AHAT 

725*  BINIT  =  BHAT 

724=  9  IFUSWANA.EQ.0)  GO  TO  10 

727=  CALL  QUANCD(XPHAT,XX,F) 

728*  GO  TO  11 

729=  10  CALL  NORMCD (X PHAT*X  X  ,KD  ,F) 

730=  11  A  1=0 
731=  A2=0 

732*  A 4=0 

733=  B1=0 

734=  B2=0 

735=  WRITE(7t52) 


73 6*  52  FORMAT(1H0,'*+  2  **0 

737*  DO  12  I  *  1,XD 

738*  WRITE<7,100)  I,XXU),XPHAT(I) 

739*  190  FQRMAT(1H0,'I,XX(I)»XPHAT(1>=  ',I5,2(F10.5» 

740*  W  *  XN(I)  /  GCPHAT(l)*(i.  -  XPHAT(I) ) ) 

741*  FF  *  F<I)*F(I) 

742*  A1  *  A1  +  W«FF 

743*  A2  *  A2  +  XX(D*W*FF 

744*  A4  *  A4  ♦  XX(I)*XX(1)*W»FF 

745*  BB  *  W*F(D*(XSAVEF(D  -  XPHAT(I) ) 

746*  B1  *  B1  ♦  BB 
747*  12  B2  *  B2  +  BB#XXO) 

748*  WR1TE<7,53) 

749*  53  FORMAK1H0,'**  3  **0 
750*  DET  -  A1*A4  -  A2*A2 
751*  WRITE<7,54) 

752*  54  FORMAT<1H0,'«*  4  **•) 

753*  DET  *  1./DET 

754*  ADELT  *  DET*(A4*B1  -  A2*B2) 

755*  BDELT  *  DET*(A1*B2  -  A2*B1) 

756*  ABSAD  *  DABSIADELT) 

757*  ABSBD  *  DABSffiDELT) 

758*  IF(  ©MAI  1  (ABSAD»ABSBD).GE.EPSI).AND.(ITERA  ©T.MAI  ID  )  00  TO  14 

759*  GO  TO  15 

760*  14  AHAT  *  AHAT  ♦  ADELT 

761*  BHAT  *  BHAT  ♦  BDELT 

762*  ITERA  -  ITERA  ♦  1 

763*C«GO  TO  ITER 

764-  GO  TO  7 

765*  15  IFf  ©MAX  1(ABSAD»ABSBD)X}EJEPS1)^HD.(ITERA.EQJ1AIIT)) WRITE  (7,16) 

766*  1MAXIT 

767*  16  F0RMAT(1H9/THE  ITERATION  PROCESS  HAS  BEEN  STOPPED.  CONVERG  ENCE  HA 
768*  IS  NOT  BEEN  ATTAINED  AFTER  ',14,'  ITERATIONS.') 

769*  A1  *  DET*A1 

770*  A2  *  DET*A2 

771*  A4  *  DET*A4 

772=  XLIXE  *  0. 

773*  DO  17  1*1  JO) 

774*  YN  *  XN(I) 

775*  YS  *  XS(I) 

776*C*+***Q1  *  DLGAMACYN  +  1 J  -  DLGAMA(YS+i  J  -  DLGAMACYN  -  YS  +1 J 
777*C»  KEHL  FIX  DLGAMA  NOT  ON  IMSL 
778*  XrEHLi*DABS<GAMMA<YN+iJ) 

779*  XXEHL2*DABS(GAMMA  (YS+i  J) 

780*  IXEHL3-DABS(GAMMA(YN-YS+U) 

781*  Q1*DL0G(XKEHL1)-DLQG0CXEHL2M)LQG<XEEHL3) 

782*  17  XLIXE  «  XLIKE  ♦  Q1  ♦  YS*DLOG(XFHAT<I))  ♦  <YN-YS)*DL0GC1.-XP  HATH)) 

783*  RETURN 
784*  END 


APPENDIX  F 


Computer  Proorams  -for  Section  VI 

The  first  two  programs,  SYMTST  and  ASYMTST,  calculate  the  test  sta¬ 
tistic  in  equations  6.3  and  6.7  respectively.  The  input  data  should  be 
on  TAPE  5  in  the  following  order: 

1. )  First  card  -  the  logit  regression  coefficients  B0  and  Bl. 

2. )  Second  card  -  title  up  10  characters. 

3. )  Third  card  -  N,  the  number  of  dose  levels. 

4. )  Fourth  through  N+3  cards  -  the  dose,  number  of  subjects,  number 

of  responses 

The  test  statistic  and  its  variance  are  both  printed  on  the  interactive 
system  display  and  written  to  TAPE  7. 

The  next  two  programs,  SYMFIT  and  ASYMFT,  uses  an  incremental  pro¬ 
cedure  to  fit  the  data  to  the  models  given  by  equations  6.1  and  6.5  re¬ 
spectively.  The  value  of  lambda  is  incremented  and  new  values  of  the 
least  squares  coefficients  are  calculated  until  the  log  likelihood  func¬ 
tion  is  maximized.  The  input  to  the  programs  should  be  on  TAPE  5  in  ex¬ 
actly  the  same  order  as  2  through  4  above,  the  first  card  is  NOT  used  by 
these  programs. 


1*  PROGRAM  SYMTST 

2*  REAL  N(20),X<20),R<20),THE<20>,IBL<2>,IBB<2,2>,ILL,TCUBE(20  ) 
3*  +»E(20) 

4*  CHARACTER  NAME#10 
5*  READ(5»*)B0fBl 
6*  R£AD(5/(A10)/)NAME 

7*  READ(5,*)NUM 
8*  DO  10  I*1|MUM 
9*  READ<5,*)X(I),N(I),R(I) 

10*  X(I)*ALOG(X(I)) 

11*  IF«(I).EQ.0)  THEM 

12*  R<I>*.5 

13*  END  IF 
14*  TEMP*B0+B  1  *X  (I) 

15*  TCUBE(I)*TEMP**3 
16*  THE<I)*l/<i.+EXP<-TEMF)> 

17*10  CONTINUE 

18*  U*0.0 

19*  ILL*0.0 

20*  DO  20  1*1,2 

21*  IBL(I>*0.0 

22*  DO  30  J*l,2 

23*  IBB(IhJ)*0.0 

24*30  CONTINUE 

25*20  CONTINUE 

26*  DO  40  1*1, NUM 

27*  FACT*N(I)*(THEU>*<1.-THE<I)» 

28*  ILL»ILL+(FACT*TCUBE(I)**2)/144 

29*  PART*(«a)~N(D#THE(I))*TCUBEU))/12 
30*  U-U+PART 

31*  IBL(l)*IBL<i)+<FACT*TCUBE<D)/i2 
32*  IBL(2>«IBL<2)+<FACT*X<I)*TCUBE(I»/i2 
33*  IBB(1»1)*IBB(1»1)+FACT 

34*  IBB(1,2>*IBB<1,2>+FACT*X<1) 

35*  IBB(2*2)*IBB(2,2)+(FACT*X  (I)**2) 

36*40  CONTINUE 

37*  DET*IBB(1 »i)*IBB(2*2)-IBB(  1 ,2)  **2 

38*C  REDEFINE  IBB  TO  BE  IBB  INVERSE 

39*  TEMP*IBB(1,1)/DET 

40*  IBB(1,1)*IBB(2,2>/DET 

41*  IBB(2,2)»TEMP 

42*  IBB(1,2)*-IBB<1,2>/DET 

43*  IBB(2,1>*IBB<1,2> 

44* C  FINISH  INVERSE  ROUTINE 

45*  Tl*04 

46*  T2*0.0 

47*  DO  50  1*1,2 

48*  T1*IBL(1)«IBB(I,1HT1 

49*  T2*IBL<I)#IBB<I,2)+T2 

50*5#  CONTINUE 

51*  VAR*ILL-(T1*IBL(1)+T2*(IBL(2))) 

52*  STDV*SQRT(VAR) 

53*  PRINT*, U^TDV 


54*  TESTV-D/STDV 
55*  PRINT*, TESTV 

56*  WRITE!?, '(A  10»3F20.10)')NAME*U,STDV,TESTV 

57*  END 


I-  PROGRAM  ASYMTS 

2*  REAL  N(20),X(20),R(20),THE(26)»IBL(2)«IBB(2,2),ILL,FACT(20) 
3*  +,E(20> 

4*  CHARACTER  NAME*10 
3*  READ(5,*)B0,B1 

6-  REAO(3»/(A10)ONAME 

7-  READ(5,*)NUM 
8*  DO  10  I-1,NUM 

9*  READ(5,*)X  (I)»N(I)»R(D 

10-  X(D«ALOG<X(D> 

II-  IF«<I).EtL0>  THEN 

12-  RCD-.5 

13-  END  IP 

14-  TEMP«B0+B1*XU> 

15-  E(I)«EXP(TEMP) 

16-  THEd)-E(I)  /  (1  .+E(I)) 

17- 10  CONTINUE 

18-  U-0.0 

19-  ILL-0.0 

20-  DO  20  1-1,2 

21-  IBL(I)-0.0 

22-  DO  30  J*l,2 

23-  IBB(W)»0.0 

24- 30  CONTINUE 

25- 20  CONTINUE 

26-  DO  40  I-14IUM 

27-  FACT(I)»THE(D+ALOG(l  .-THE(I» 

28-  ILL»ILL-KFACT(I)**2)*N(I)/E<I) 

29-  U»U+(FACTa)/THE(I))»<R(D-N(I)#THEa)) 

30-  IBL(l)-IBL(l)*PACT(I)*Na)*(l.-THE(I)) 

31*  IBL(2)-IBL<2)+FACT(I)*N(I)*X(I)*d.~THE(I)) 

32-  IBBd»l)-IBBd,l)+N(D*THE(I)*d.-THE(I)) 

33-  IBBd  ,2)*  IBBd  ,2)+N(D*X  (I)*THE(I)*d  .-THE(I)) 

34-  IBB(2,2>«IBB<2,2)*Na)*<X(I)**2)*THEa)*d.-THE<I)) 

3$iu  rnuTunre 

36-  DET-IBBd  ,l)*IBB(2»2)-IBBd  ,2)**2 

37- C  REDEFINE  IBB  TO  BE  IBB  INVERSE 

38-  TEMP»IBB(1,1)/DET 

39-  IB8d,l)»IBBC2,2)/DET 

40-  IBB(2»2)«TEMP 

50-  IBBd  ,2)— IBBd  »2)/DET 

51-  XBB(2»1)«IBB(1»2) 

52- C  FINISH  INVERSE  ROUTINE 

53-  T1-0J 

54-  T2-0J 

55-  DO  50  1-1*2 

56-  Ti-IBL(I)*IBBa,l>*Ti 

57-  T2-IBL(I)*IBBd,2)+T2 

58- 50  CONTINUE 

59-  VAR-ILL-(Tl*IBLd)+T2*(IBL(2))) 

60-  STDV*SQRT(VAR) 

61-  PRINT*, U3TDV 

62-  TESTV-U/STDV 


1*  PROGRAM  SYMFIT 

2*  IMPLICIT  DOUBLE  PRECISION(A-H,0-Z) 

3*  COMMON/SET1/X(20),XP<20>JAU(20),XLAMDA,AHAT,BHAT,N(2&) 
4*  COMMON/SET2/NUM.1XSW 
5*  REAL  ISTAT>IDF»IQUE 

6*  DOUBLE  PRECISION  MLEHAT,R(2«),PFINC28),THE(20) 

7*  CHARACTER  NAME«20 
8*  1XSW*1 
9*  CHECK*  10000. 

18*  CHXLII*-999999999. 

11*  PRINT*/ NAME?' 

12*  READ  (*/(A20)  ONAME 

13*  R£AD(5,*)NUM 
14*  IDF*NUM-2. 

IS*  DO  10  1N*1,NUM 

16*  READ(5i*)X(IN>tN(IN)rR(IN) 

17*  X  (IN)*DLOG(X  (1N» 

18*  IF(R(IN)JBQ.  0JRUN}*.5 

19*  XP<IN)*R(IN)/N(IN) 

20*10  CONTINUE 

21*  INCLAM*0 

22*  DO  20  J*l>200 

23*  X  LAMDA*DBLE(INCLAM)  / 1 00. 

24*  INCLAM*INCLAM+1 

25*  XLLF*0.0 

26*  CALLLSE 

27*  STAT*0 

28*  XLIKE*0.0 

29*  DO  30  I*1»NUM 

30*  FACT*J5*ILAMDA*(AHAT+BHAT*XU)) 

31*  IF®ABS(FACT)J.T.  1.  -AND.  FACT  .NE.  0.)  THEN 

32*  T1»(1.+FACT)**<1/XLAMDA) 

33*  T2*(1.'FACT)**(1/XLAMDA) 

34*  THE(I)*T1/(T1+T2) 

35*  ELSE  IFIFACT.EQ.  0.0)  THEN 

36*  THE(I)*i/<l+EXP(-AHAT-BHAT*I<I>» 

37*  ELSE  IFtFACT-LE.  -1J  THEN 

38*  THE(I)*.9999999999 

39*  ELSE 

40*  THE(I)*.0000000001 

41*  ENDIF 

42*  EXPECT*N<I)*THE(D 

43*  DIFFSQ*(R(I)-EXPECT)#*2 

44*  TEST*DIFFSQ/(EXPECT*(1.-THE(I») 

45*  STAT*STAT+TEST 

46*  XLIKE»XLIIE+Ra)«DLOG(THE«))+(N(I)-R(I))*DLOG<l.-THEa)) 

47*30  CONTINUE 

48*  IF(STATXT.CHECDTHEN 

49*  DO  40  XX*1»NUM 

50*  PFINOCX)*THE(XX) 

51*40  CONTINUE 
52*  CHECK-STAT 

53*  CHXXIX-XLIEE 


54*  MLEHAT-XLAMDA 

55*  B4-AHAT 

54*  Bl-BHAT 

57*  END  IF 
58*24  CONTINUE 

59*  WRITE(7»'(F4«4»2F24.15)')MLEHAT»B4'B1 

44*  WRITE<7,«CHILIX 

41*  MSITE(7,'(A24)')NAME 

42*  WRITE!? »*) THERE  ABE  *,NUM,*  DATA  POINTS.' 

43-  HRITE(7,«) 

44*  WRITE  (7,*) 

45*  WRITEf?,*) 

WRITE  (7>*)'  MTBF  LAUNCHES  ABORTS  PREDABT  FROB 

47*  DO  54  I*1»NUM 

48-  EXPECT*N(D«PFIN(D 

49*  DIFFSQ«<R(I)-EXPECT>**2 

74*  TEST*DIFFSQ/<£XPECT*<i.-PFIN(D» 

WRITE(7»'(3F7.0,7X  »F?.5»5X  »F4.5»3X  »F7 .4)')X  (I),REAL(N(I) 

72*  +,£XPECT,PFIN(I),TEST 

73*54  CONTINUE 

74*  ISTAT*REAL(CHECX> 

75*  CALL  MDCH(ISTATflDF»IQUE»IER) 

74*  PCHI*1.-IQUE 
77*  WRITE(7,*) 

78*  WRITEC7,*)'  TEST  STATISTIC  CHI-SG  TAIL  PRDB' 

79*  WRITE(7t'(3X«F14.6»15X»F7.5)')CHECKrPCHI 
84*  WRITE(7,«) 

81*  END 

82*  SUBROUTINE  LSE 

83-  IMPLICIT  DOUBLE  PRECISION(A-H,0-Z) 

84-  COMMON/SET1  /X  (28)»XP(29)rTAU(24)»XLAMDA  »AHATiBHAT»N(24) 
85*  COMMON/SET2/NUK.IXSW 

84*  DOUBLE  PRECISION  W<24> 

87*  DO  1  I*1*MUM 

88*  IFGCLAMDA.EQ.  0.0)  THEN 

89*  TEMP*XP<I)/(1.-XF<I)> 

94-  TAU(I)*DLOG(TEMP) 

91*  ELSE 

92-  TEMP1*<1.-XP<I))*#XLAMDA 

93*  TEMP2*XP<I)**XLAMDA 

94*  XNUM*2#<TEMP2-TEMP1> 

95-  XDENOM«XLAMDA*(TEMP2+TEMPl> 

94-  TAU(l)*XNUM/XDENOM 

97*  END  IF 
98*1  CONTINUE 
99*  IFIIXSW  .EQ.  0J  GOTO  4 
144-  SUMX*4. 

141*  SUMW-4.4 

142*C  SET  UNW-4  IF  UNWEIGHTED  LS  IS  WANTED  SET  TO 
143*C  ANYTHING  ELSE  FOR  WEIGHTED  LS 
144*  UNW*1. 


105*  ITX*0. 

106*  DO  5  1=1, NUM 
107*  WU>*DBLE(Na)>*XP<I)*{l.-XPa» 
108=  IF(UNW.EQ.0.)THEN 

109=  W(D*i. 

110=  SUMW=SUMW+W(I> 

111=  ELSE 

112=  SUMW*SUMW+W(I> 

113*  END  IF 

114=  SUMX*SUMX+X(I)*W(I) 

115*  XTX«XTX+CC<I>**2)*W<I> 

114*5  CONTINUE 

117*  SSX*XTX-<SUMX»2>/SUMW 

118=  X  BAR3  SUM!  /  SUMW 

119*  IXSH-0 

120=4  CONTINUE 

121*  ITY*0. 

122=  SUMT*0. 

123=  DO  6  1*1, NUM 

124=  XTY*XTY+X(I)*TAU«)*W<I> 

125*  SUMY  =SUMY +TAU(I)*W(I) 

124=6  CONTINUE 

127*  Y  BAR*  SUMY /SUMW 

128=  SSP=XTY-(SUMX*SUMY)/SUMW 

129*  BHAT=SSP/SSX 

130=  AHAT=YBAH-BHAT*IBAR 

131*  RETURN 

132=  END 


135*  PROGRAM  ASYMFT 

136=  IMPLICIT  DOUBLE  PRECISION(A-H,D-Z) 


137*  COMMON/SETl/X(20),XP(20),TAU(20),XLAMDA,AHAT,BHATtN(20> 
138=  COMMON/SET2/NUM,IXSW 
139=  REAL  1STAT,IDF,IQUE 

148*  DOUBLE  PRECISION  MLEHAT,R(20),PFIN(20>,THE(20> 

141=  CHARACTER  NAM£*60 

142=  IXSW=i 

143-  CHECK*  10888. 

144*  CHKLIK*-9999999999 

145*  PRINT*,'  NAME?' 

146=  READ(*»'(A60)')NAME 

147=  READ(5,*)NUM 
148=  IDF=NUM-2. 

149-  DO  10  IN=1,NUM 

150=  READ(5,»)X(IN),N(IN),R(IN) 

151*  X  <IN)*DLOG(X  (IN)) 

152=  IF(R(IN).EQ.  0JR<IN>=.5 

153*  XP(IN)*R(IN>/N<IN) 

154=  IF(XP<IN).Ea.l  )XP(IN)=.9999 

155=10  CONTINUE 

156=  INCLAM—500 

157*  DO  20  J=  1,1000 

158=  XLAMDA-DBLEdNCLAM)/ 1. 

159*  INCLAM*  INCLAM+ 1 

168=  XLLF*8.8 
161-  CALL  LSE 
162=  STAT=8 
163-  X  LIKE-0.0 

164*  DO  30  I-1.NUM 

165-  FACT=XLAMDA*EXP(AHAT+BHAT*X(I)) 

166=  IF(FACTX>T.  -1.  .AND.  FACT  .NE.  0.)  THEN 
167*  T1=(1.+FACT)**(-1./XLAMDA) 

168*  THE(I)*1.-T1 

169*  ELSE  IF<FACT.EQ.  0.0)  THEN 

170=  THEd)*i.-EXP(-EXPCAHAT+BHAT*X(D» 

171-  ELSE  IFCFACT.LE.  -1.)  THEN 

172*  THE(I)=  .9999999999 

173*  END  IF 
174=  EXPECT=Nd)*THEd) 

175*  DIFFSQ*  (R  (I)-EX  PECT)  **2 

176=  TEST’D  IFFSQ/(EXPECT*(1.-THE(I)W 

177=  STAT*STAT+TEST 

178=  XLIKE*XLIKE+RU)*DLOG(THE(I))+QI(I)-R«))*DLOG(l.-THE(I)) 

179=30  CONTINUE 

180*  IFCXLIKE.GT.CHKLI10THEN 

181*  DO  48  KK*1,NUM 

182*  PFIN(KK)*THE(KK) 

183-40  CONTINUE 
184*  CHECK-STAT 

185-  CHKLIK*XLIKE 

186-  MLEHAT*XLAMDA 


187*  B0*AHAT 

188*  B1*BHAT 

189*  END  IF 
198*20  CONTINUE 

191*  WRITE(7,'<F10.5,2F20.15)')MLEHAT,B0,Bt 

192*  WR1TE(7  ,*)CHKLI  K 

193*  PRINT*, MLEHAT 

194*  NR1TEC7  ,'(A60)QNAME 

195*  NRITE(7,*)'THERE  ARE  ',NUM,'  DATA  POINTS.' 

196*  WRITEC7,*) 

197*  NRITEC7,*) 

198*  NRITEC7,*) 

^*  WRITEf?,*)'  MTBF  LAUNCHES  ABORTS  PRED  ABT  PROB 

208*  DO  50  I*1,NUM 
201*  EXPECT*N(I)«PFIN(I) 

202*  DIFFSQ*(R(1)~EZPECT)*»2 

203*  TEST*DIFFSQ/(EXPECT*(1.-PFIN(I))) 


m 

205* 

206*50 

207* 

208* 

209* 

210* 

211* 

212* 

213* 

214* 

215* 

216* 

217* 

218* 

219* 

220* 

221* 

222* 

223* 

224* 

225* 

226* 

227* 

228* 

229*1 

230* 

231* 

232* 


NRITE(7  ,'(3F7 .0 ,71  ,F9 .5 ,5X  ,F6 .5  »3X  »F7 .4)  T),REAL(N(I) 

+, EXPECT, PFINd), TEST 
CONTINUE 
ISTAT»REAL(CHECX) 

CALL  MDCH(ISTAT,IDF,1QUE,IER) 

PCHI*1.-IQUE 

WRITEC7,*) 

WRITEC7,#)'  TEST  STATISTIC  CHI-SQ  TAIL  PROB' 

WRITE(7,'(3X,F10A,i5X,F7^)OCHECK,PCHI 

WRITEC7,*) 

END 

SUBROUTINE  LSE 

IMPLICIT  DOUBLE  PRECISIQN(A-H,0~Z) 

COMMON/SET1/X(20),XP(20),TAU(20),XLAMDAAHAT,BHAT,N(20) 

COMMON/SET2/NUM,IXSW 

DOUBLE  PRECISION  H(20) 

DO  1  1*1, NUM 
IFCXLAMDA.EQ.  0.0)  THEN 
TEMP*l./(l-XPa» 

TAU(I)*DLOG(DLOG(TEMP» 

ELSE 

TEMPl*(i.-XP<I))**<-XLAMDA> 

TEMP2*(TEMP1-1J/XLAMDA 

TAU«)*DLOG(TEMP2) 

END  IP 
CONTINUE 

IFUXSH  .EQ.  0J  GOTO  4 
SUMX*0. 

SUMN*04 


233-C  SET  UNW*0  IF  YOU  WANT  UNWEGHTED  LS,IF  NEED  WEIGHTED 


234*C  SET  UNW  TO  ANYTHING  BUT  1  (ONE) 


235*  UNH*0. 

236*  XTX*0. 

237*  DO  5  1*1, NUM 


238=  W<I)=DBL£(N(I))»XP(1)*(1.-XP(I)> 

239=  IP(UNW.EQ.8JTHEN 

248=  W<I)*1. 

241=  SUMW=SUMW+W<I) 

242*  ELSE 

243=  SUMW=SUMW+W(I) 

244=  END  IF 

245*  SUMX=SUMX+X(I)*W<I) 

246*  XTX*XTX+<X<I)**2>»W<I) 

247=5  CONTINUE 

248=  SSX=XTX-<SUMX#*2)/SUMW 

249*  IBAR-SUMX/SUMW 

258*  IXSW*8 

251=4  CONTINUE 

252*  ITY=8. 

253=  SUMY =8. 

254=  DO  6  1=1, NUM 

255=  XTY=XTY+X<I)*TAUU)#W(I) 

256*  SUMY=SUMY+TAUU)*W(I) 

257=6  CONTINUE 
258*  YBAR=SUMY/SUMW 
259*  3SP=XTY-<SUMX*SUMY)/SUMW 
268*  BHAT*SSP/SSX 
261=  AHAT=YBAR-BHAT*XBAR 
262*  RETURN 
263* 
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available  techniques.  An  alternative  to  performing  numerous,  and  some¬ 
times  redundant,  simulations  is  to  use  these  techniques  whenever  pos¬ 
sible. 

Data  from  the  Avionics  Evaluation  Program  (AEP>  were  used  as  the 
basis  for  estimating  the  probability  of  aircraft  abort,  based  on  the 
mean- time-between -failure  (MTBF)  of  various  equipment  items,  using  four 
quantal  assay  techniques.  The  fits  obtained  from  these  models  were  com¬ 
pared  to  the  more  popular  probit  and  logit  results  previously  obtained 
by  Dr.  David  Barr. 
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